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Abstract 

During the last decade, many exciting phenomena have been experimentally observed and theoretically 
predicted for ultracold atoms in optical lattices. This paper reviews these rapid developments concentrating 
mainly on the theory. Different types of the bosonic systems in homogeneous lattices of different dimensions 
as well as in the presence of harmonic traps are considered. An overview of the theoretical methods used 
for these investigations as well as of the obtained results is given. Available experimental techniques are 
presented and discussed in connection with theoretical considerations. Eigenstates of the interacting bosons 
in homogeneous lattices and in the presence of harmonic confinement are analyzed. Their knowledge is 
essential for understanding of quantum phase transitions at zero and finite temperature. 
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1. Introduction 


Ultracold atoms in optical lattices opened a new era in the study of quantum many-body phenomena. In 
contrast to other condensed-matter systems, they provide a unique opportunity of control. Using interference 
of laser beams propagating in different directions, one can create various types of periodic potentials with 
the amplitude proportional to the laser intensity mm that are free of defects and dissipative channels. 

Optical lattices provide an efficient tool to control the system dimensionality. Apart from the three- 
dimensional geometry [ 3 ] , it is possible to reach very strong spatial confinement in certain spatial directions 
and reduce system dimensionality creating single or periodically arranged linear HD and planar |7HI3 
structures. Rapidly moving laser beams allow to arrange ring lattices and two-dimensional periodic struc¬ 
tures with point-like defects m- Disorder with known statistical properties and tunable parameters can be 
also introduced into the system [HHn] either by optical means using incommensurate optical lattices |18| 
and laser speckles m or through the interaction with other atomic species localized at random positions |20| . 

Most of experiments in optical lattices are performed with alkali-metal atoms, mainly with Rb (see, 
e.g., PHSl[3mHlS] ) and also with Li j25J|27j) Na [231152], K [20], Cs |B1[T21ISI|- In recent years this list was 
extended and includes also Yb [ 55 ], which belongs to the atoms of the alkaline-earth metals, and Cr [ 55 ] 
which is a transition metal. The choice of the atoms is mainly determined by the fact that their electronic 
transitions lie in a convenient spectral range allowing efficient manipulation by an optical laser. The atoms 
can be trapped either in a single state or in a manifold of the electronic ground states and cooled below the 
temperature of quantum degeneracy. If the total number of electrons, protons and neutrons which constitute 
the atoms is even, the latters are bosons, otherwise they are fermions. In this review, we shall consider only 
bosons. 

Two-body interactions of these atoms, except Cr, are of short range and the effective strength can be 
controlled by the intensity of the laser creating the optical lattice or by Feshbach resonances |34fl36| . The 
latter is accompanied by the formation of molecules that are converted back to atoms. 

In experiments, ultracold atomic system can be controlled on a macroscopic as well as microscopic levels. 
Macroscopic measurements are based on the time-of-fiight imaging [3] El] and Bragg spectroscopy [58MT| 
which provide information about the energy spectrum and the state of the system in the momentum space. 
More recently, new techniques have been developed to perform in situ measurements on a microscopic 
level m with the spatial resolution of the order of one lattice period or even less. Tremendous progress 
has been also achieved in the single-site and single-atom addressability [EHinjiliiiiii which is important for 
applications in quantum technology. 

In deep periodic potentials, atoms can move from one potential well (lattice site) to the next one by 
quantum tunneling which gives rise to the discrete lattice models. In the case of one-component spinless 
bosons, the lattice system is described by the Bose-Hubbard model. This model was originally introduced in 
a rather heuristic manner in order to describe the differences in the ground state and low-energy excitations 
of interacting bosons in a homogeneous space under density variations and the associated solid-superfiuid 
transition in '^He |431l^ . Later it was derived in the context of the solid state physics @6] and motivated by 
experiments on ^He absorbed in porous media or Cooper pairs in granular media m- The presence of spin 
degrees of freedom and Feshbach resonances lead to extensions of the standard Bose-Hubbard model. In the 
case of cold atoms, the parameters of the corresponding lattice model can be derived from first principles 
which allows direct comparison of the theoretical predictions with experimental data. 

A remarkable feature of the Bose-Hubbard model is that it reveals a quantum phase transition from 
the superfiuid to the Mott insulator [UliH] that results from the competition between the kinetic energy 
and on-site interaction. It is characterized by a natural order parameter - the superfiuid fraction. In the 
case of spinless bosons, it is a second-order transition. In the superfiuid phase, the spectrum of excitations 
has no gap and the particle-number statistics is described by a broad Poisson-like distribution. In two and 
three dimensions, the one-body density matrix (two-point correlation function of the first order) shows the 
off-diagonal long-range order and decays as a power law in one dimension. In the Mott-insulator phase, there 
is a finite gap in the excitation spectrum, particle-number fluctuations are suppressed, and one-body density 
matrix decays exponentially in all dimensions. In low dimensions, Mott insulator possesses nonlocal string 
order ^21 [ 50 ] . Superfiuid-Mott-insulator transition in optical lattices has been experimentally observed first 
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in three dimensions [S] and then in one and two dimensions. The presence of spin degrees of freedom and 
Feshbach resonances leads to qualitatively new features. 

In spite of a big progress in theory, complete description of the interacting quantum systems is still a 
challenge. Though at first glance seemingly very simple, even the standard Bose-Hubbard model is not 
analytically solvable in general. Exact analytical solutions are known only in very special situations, as the 
case of vanishing or infinitely strong interaction. Approximate analytical results are obtained by systematic 
expansions in powers of small parameters. However, they have always their limitations. For instance, strong¬ 
coupling expansion |5HI54| is valid in arbitrary dimensions and for arbitrary filling factors but limited to 
small tunneling rates. In addition, it can be easily implemented only for fillings close to commensurate due 
to the degeneracy of the superfluid state. The expansion in powers of the inverse filling factor |55[ 155] is 
valid in arbitrary dimensions and for arbitrary tunneling rates but cannot be applied if the filling factor is 
of the order of one. 

Parallel to the analytical studies, different exact numerical methods have been developed for the analysis 
of the Bose-Hubbard model. The most straightforward and easiest to implement is exact diagonalization |57F 
However, due to exponential growth of the Hilbert space with the system size, the method can be used 
if the number of particles N and the number of lattices sites L are rather small. The largest system of 
bosons reported in the literature was L = iV = 18 [SS], although with the restriction that no more than four 
particles can occupy one lattice site. 

More sophisticated deterministic method is the density-matrix renormalization group uniiii] which is 
based on the fact that usually the state of the system occupies only a small subspace of the exponentially 
large Hilbert space. This approach is quite successful in one dimension and allows to treat larger systems {N 
and L of the order of 1000 m) but it fails in higher dimension, where quantum Monte Carlo methods CSJEll 
become more efficient. By stochastic sampling they allow to treat stationary states in realistic experimental 
situations of TV = 3 x 10® bosons in a three-dimensional lattice |24| . 

Mean-field theory plays an important role in the studies of the Bose-Hubbard model. It is based on 
the Gutzwiller ansatz m which takes into account local fluctuations but neglects quantum correlations of 
different lattice sites. This approach is exact in infinite dimensions and provides a useful insight into the 
physics in large finite dimensions but it fails for low-dimensional systems. Attempts to correct the mean-field 
theory incorporating distance-dependent quantum correlations were undertaken by several authors. These 
include random phase approximation US El], cluster mean field m, method of effective potential 1291180], 
dynamical mean-field theory (DMFT). The latter is probably the most successful approach but computa¬ 
tionally quite demanding. It was originally developed for fermions [ST] and recently for bosons |82H86| . It 
can be derived as an expansion in powers of the inverse coordination number |86| and reduces to the solution 
of an impurity problem on a single site or a cluster of sites. This is a difficult computational problem which 
requires application of exact methods like exact diagonalization, DMRG, QMG. 

In recent years, excellent reviews on cold atoms in optical lattices were published |87ti55] . However, many 
important aspects were not properly discussed and the field continues to grow. The plan of this review is 
the following. In sectionwe discuss the basic mechanism for the creation of external potentials for neutral 
polarizable atoms by optical laser fields and consider eigenstates of single atoms in periodic potentials of 
different types. This serves as a preliminary step for the derivation of the Bose-Hubbard models of various 
types. In section we derive the Bose-Hubbard Hamiltonian for the simplest case of spinless bosons and 
discuss its symmetries. In section we provide definitions of basic physical quantities that are used in 
the theory of low-temperature phenomena in a lattice. Section provides theoretical background of the 
main experimental techniques for cold atoms in optical lattices. In section we present exact results for 
bosonic many-body systems in the simplest special cases that allow analytical treatments. Section gives 
an overview of the perturbative results for the ground state and lowest excited states in the regime of strong 
interactions. Section is devoted to the criticality of the Bose-Hubbard model and the quantum phase 
transition from the superfluid to the Mott insulator. In section we present exact numerical results for 
macroscopic and microscopic quantities across the quantum critical point. In section |10[ we review the 
mean-field theory based on the Gutzwiller approximation. In section 11 we consider a system of lattice 


bosons near a Feshbach resonance. Section [T^ deals with physics of spin-1 bosons. 
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hyperfine splitting 


fine splitting 



Figure 1: Scheme of the electronic levels for atoms with nuclear spin 7 = 3/2 (not to scale). 


2. Single atom in a periodic potential 

Alkali atoms consist of a spherically symmetric atomic residue and one outermost electron (spin S = 1/2) 
in the state with a principal quantum number n, that maybe different for different atoms, and orbital angular 
momentum L. The full scheme of the relevant electronic levels for atoms with nuclear spin I = 3/2 as in the 
case of ^Li, ^^Na, and ®^Rb, is shown in Fig. The ground state is an S-state (L = 0). Spin-orbit 

coupling leads to the fine splitting of the first excited level (P-state with L = 1) into two states separated by 
the energy Aps- The states are distinguished by the values of the electronic angular momentum J = L + ^ 
and form the D-line doublet ^Si /2 —^Pi /2 (Dl), ^Si /2 —^P 3/2 (D2)g The coupling of the electronic 
spin to the nuclear spin then leads to the hyperfine splitting of both ground and excited states with the 
energies fiAupS) ^^hfsi ^^hfs- additional coupling provides hyperfine levels with the total angular 
momentum (hyperfine spin) F = J -|- I which are manifolds of 2F + 1 degenerate states characterized by 
the magnetic quantum numbers mp = 0, ±1,..., FF. The energies of the hyperfine splittings are five orders 
of magnitude smaller than for the fine splitting. 

Laser field acting on the atom causes different transitions between electronic levels, which are determined 
by the frequency and polarization of the laser wave. The transitions from ^Si /2 to ^Pi /2 and ^P 3/2 are electric 
dipole transitions. They are allowed if the selection rules Amp = 0 for linear polarization or Amp = ±1 
for circular polarization are fulfilled [M]. If the detuning of the laser frequency is much larger than the 
spontaneous emission rate, one can adiabatically eliminate all the excited states in the spectrum of atoms 
denoted by F' and F" in Fig.[^ This leads to the effective potential acting only on the ground state sublevels 


^For given L and S, J takes the values in the range J = \L — S\,L + S. 

^Here we use a standard spectroscopic notation for the states with the principal quantum number n, orbital 

angular momentum L = S, P,..., spin S, and electronic angular momentum J. 

®In analogy to the total angular momentum J, F takes the values in the range F = | J — 7|,..., J -|- 7. 
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labeled by F (see, e.g., |951 [55] and references therein) 

V- (EW ■ d7a)* (E(x) • ^ 

hiuJL-UJ,) 

where E(x) is the electric field strength of the laser field, is the dipole matrix element between the 
gronnd state snblevel a and the excited state snblevel 7 of energy This allows to create controlled 

potentials for nentral polarizable atoms which can be of completely different types ranging from random and 
quasi-random to perfectly periodic. In this review, we will be dealing with the potentials of the latter type 
known under the name optical lattices. The overview of the geometries of the optical lattices was given in 
Refs. [TJ [21 EZl IHE] . Here we consider mainly hypercubic lattices but focus more on the effects coming from 
the interference of the excited electronic levels which have various manifestations depending on the laser 
frequency as well as polarization. 


2.1. One-dimensional lattice in the case of large detuning 

We consider a pair of counterpropagating laser beams with the wavevectors Icl and —Icl along the Xi- 
direction. If the detuning is much larger than the hyperfine splitting of the electronic levels, this laser 
configuration does not lead to any coupling of the internal ground states. It creates a one-dimensional 
periodic potential which is the same for all ground-state sublevels and has the form 

Vl{x) = VoCOs‘^ , ( 2 ) 

where a = Tr/fen = Al/ 2 is the period (lattice constant). If the two laser beams intersect at an angle ip < tt, 
one can create a one-dimensional lattice with a larger period given by a = Al/ {2 /2)). Using this 

technique, optical lattices with a up to 80 ^m were demonstrated in experiments with ®^Rb in the field of 
Ti:Sa laser emitting at the wavelength Al = 820 nm |99| . 


2.1.1. Bloch bands 

We suppose that the system consists of L potential wells and impose periodic boundary 
the wavefunction of the atom ijjlx La) = iIj{x) which satisfies the Schrodinger equation 




+ Vi,{x) 


tpix) = E'ip{x) , 


X e 


La La 

T’ T 


According to the Bloch theorem, the solution has the following form 

'tp{x) = 'tpb{x; k) = Ub{x\ , E = Eb{k) 


conditions on 


(3) 

(4) 


where Ubfx] k) is a periodic function of x with the period a and b is the band index. The wavenumber 
k = kq = 2TTq/{La) takes in general discrete values determined by the integer q which is defined up to 
modulo L. If we do not want to care about the differences between even and odd L, we can assume that 
q = 0,..., L — 1. In this case, fc = 0,..., 2tt{L — 1) /{La). However, usually the first Brillouin zone (IBZ) is 
defined as k € [—tt/o, tt/u] and we will also follow this convention. In the limit of infinite lattice {L —)■ 00 ), 
k becomes a continuous variable. 

The solution of the eigenvalue problem can be expressed in terms of Mathieu functions. Despite they are 
rather well studied in the mathematical literature |100fll02] , exact results can be obtained only numerically. 
One can use Mathematica (see, e.g.. Ref. |103| for notes on that) but in order to have full flexibility it is 
better to write an own program, for instance, in C/CH—h or Fortran. With this purpose in mind we use the 
Fourier series expansion 


Ub{x; k) 


^ 00 

Cb„(fc)exp (f27rn-) , 
\/a V a/ 


(5) 
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where the coefficients Cb„ are the solutions of the eigenvalue problem 

OO 

'Hnn'{k)cbn'{k) = Eb{k)cbn{k) , U = -OO, . . . , OO , 

n' — — oo 


'Hnn'{k) 


ka \ Vq 
En[-+2n) +f 


Vo 

^nn' ~\~ — l E ^^n,n' + l) 


( 6 ) 


where E^, = h^k"^/{2M) is the recoil energy. They satisfy the orthonormality condition 


OO 

X] Cb^n{k)cb^u{k) = 6b,b2 ■ (7) 

n——oo 


The solutions of Eq. (|^ are periodic functions of k with the period 27r/a. This is because the shift of the 
wavenumber k —>■ fc + 27r/o can be compensated by the corresponding shift of the index n —)■ n — 1. Using 
this property one can show that ipbix] kq) and ipbix] k^) are orthogonal, unless (g' — q)/L is an integer. 

Since matrix 'H(fc) is real and symmetric, the coefficients Cbnik) can be chosen to be real which provides 
unique solutions. In addition, it is tridiagonal and the eigenvalue problem can be solved numerically using 
efficient algorithms m- In addition, the off-diagonal terms in H(fc) coincide with the lowest-order approx¬ 
imation of the second derivative via a finite difference and the diagonal terms are the same as a discrete 
harmonic potential. Therefore, one can expect that the coefficients Cbn decrease exponentially with |n| and 
the infinite-dimensional matrix 'H(fc) can be safely truncated to a moderate finite dimension. The fact that 
the system is finite leads only to the discretization of k but does not change the values of Cbnik) and Eb{k). 

Energy spectrum Eb{k) within the first Brillouin zone, is shown in Fig. At each value of k the 
spectrum is discrete and all the eigenvalues are distinct The functions Eb{k) take their extremal 

values at fc = 0, En ja. A: = 0 is a minimum for even b and maximum, if b is odd. With the increase of the 
amplitude of the periodic potential Vq, the energy bands Eb{k) become more flat and the gaps between the 
bands grow. In the limit Vq ^ Eji, the width of the bands is given by the asymptotic expression m 


\Eb{Tr/a)-Ebm 2^'^+^ 


1 - 


Er 

6&2 + Ub + 7 
16 ^ 


bl^/n 
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Figure 3; (color online) Eigenstates of Eq. § for the lowest two bands with 6 = 0 (a) and 6=1 (b). In both panels, 
Vq = 10-Er and ka/ir = 0 (red), +1 (green), —1 (blue). The lines are guide to the eye. 


The coefficients Cb„(/c) which represent the eigenstates in Eq. ([^ are shown in Fig. for the lowest two 
bands (5 = 0,1). For 6 = 0 and k = 0, Chn is an even fnnction of n. However, if we move towards the 
edges of the Brillouin zone, this symmetry is destroyed. Dne to the analogy of Eq. to the Schrodinger 
equation for the harmonic oscillator mentioned above, one would expect that con should be positive. The 
fact that con take negative values for odd n’s is simply because we have chosen Vb > 0. In the opposite case 
(Vb < 0), con are indeed always positive. In the next energy band (6 = 1), C{,„(0) is an odd function of n. 
cinik) becomes symmetric with respect to n = ±0.5 for ka/n = =Fl. Similar features can be observed in the 
higher energy bands. 


2.1.2. Wannier functions 

Bloch functions ipb{x; k) are extended over the whole lattice for any 6 and k. An alternative basis suitable 
for the description of single particles at individual lattice sites is provided by Wannier functions defined via 
the Fourier transform |106| 

Wbeix) = Wb{x -xe) = ^ ^ , (9) 

feGlBZ 


where xg = xq + a£, with £ being an integer, are the minima of the periodic potential, xq = aj^. if Vb in 
Eq. ([^ is positive but xq = 0 for negative Vb- The summation in Eq. ([^ is over the values of k within the 
first Brillouin zone. The functions satisfy the orthonormality condition |107| 



Wb,/i,{x)Wb^e^{x)dx 




( 10 ) 


and form a complete set m- They possess the symmetry Wb{—x) = {—l^Wb^x). In finite lattices, Wb{x) 
are periodic functions: Wb{x + La) = Wb{x). In the limit of infinite lattice, the sum in Eq. ([^ can be 
replaced by the integral: 


1 

L 


E - 

fcelBZ 


a 

27T 



dk . 


As it was proven in Ref. uns] for a general case of separated energy bands in one dimension, the Wannier 
functions are uniquely defined by their symmetry properties and asymptotic behavior at large distances (see 
the discussion below) that guaranties their minimal width. In the following we shall consider the Wannier 
functions for the lowest Bloch band Wq{x). 

In the limit of vanishing potential (Vb —t 0), the eigenvalue problem ([^ has a very simple analytical 
solution which leads to the following result for an infinite lattice m-- 


1 sin(7ra:/a) 
\/a TTx/a 
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Wo(cr) 


( 11 ) 




























xja xja 

Figure 4: Wannier function in the first Bloch band for Vq = 5 -Er. Solid line is exact result and dashed line is a Gaussian 
approximation (b) is the same as (a) but with a logarithmic scale. 





Figure 5: Overlap of the Wannier function for the lowest Bloch band with the Gaussian approximation | |15[ |. (V’holVFo) 
vanishes for Vq = 0 and grows rapidly with the increase of Vq . 


This function oscillates with the amplitude decreasing with the distance x, and this type of behavior is 
typical for the Wannier functions (see Fig. |^. We would like to stress that Wq{x) is not a ground-state 
eigenfunction of any Hamiltonian and the nodes appear to be necessary in order to satisfy the orthogonality 
condition (lOl. 

At finite Vq, Ijx decay of the envelope of the function Wq{x) is preserved only for |a:| <C Xc- For |a;| ^ Xc, 
the asymptotics of the envelope acquires a different form |10HL 1111] : 

:|“^/^exp (-/io|ai|) . 


Woix) 


( 12 ) 


The crossover distance Xc, which is infinite for Vq = 0, becomes finite for nonvanishing Vq and decreases with 
Vq. /iQ is a constant which vanishes in the limit Vq —)■ 0 and grows with the lattice depth Vq. For shallow 
and deep lattices it can be calculated analytically uni and the result reads |112| 


hoa _ f Vq/Er for Vq Er , 

TT “1 forl/o»i^R. 


(13) 


In the 
frequency 


case of a deep optical lattice, each lattice site can be described by a harmonic potential with the 


mi 


^ho — 


^Er ^ 

h V Er 


(14) 


Then the solution Q of the Schrodinger equation (|^ can be approximated by the eigenfunctions of the 
harmonic oscillator. This leads to the Gaussian approximation for the lowest-band Wannier function 


Wq{x) « ■i/'ho(a:) 



(15) 
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of the width 


^ho — 


Mwho 




1/4 


(16) 


We compare it with exact Wo{x) in Fig. One can see that the Gaussian approximation overestimates 
the maximum height and fails to reproduce the detailed structure of the Wannier functions. In addition, it 


violates the orthogonality condition (101 at any finite Vq. The Gaussian approximation becomes exact only 


in the limit of infinitely large Vb, when (|15| takes the form of the b-function. Nevertheless, the overlap with 
the exact Wannier function 


(V'holITo) = / 'iIj^^{x)Wo{x) dx 


(17) 


is close to one even if the amplitude Vb of the periodic potential is of the order of few recoil energies (see 
Fig. i and Ref. my Due to its simplicity, the Gaussian approximation is often used in order to obtain 
analytical estimations of the parameters of the Bose-Hubbard model |1121111414122] . 

2.1.3. Tunneling matrix 

In the basis of the Wannier functions the single-particle Hamiltonian can be represented in the form of 
the tunneling matrix with the matrix elements defined as 


J-La/2 




(a;) rfx . (18) 

Using the definition of the Wannier functions and the Bloch theorem, one can show that the matrix elements 


depend on the distance s = |^i — £ 2 ] and do not vanish only for &2 = = b- Eq. (18) can also be rewritten 

in the form 

Js = ^ ^ Eb{k) exp (ikas) , (19) 

fceiBZ 

where we observe that is a Fourier transform of the band structure. Jq is the average energy of the 
band. 

Numerical calculations show that for s > 0 the sign of alternates with the distance s and with the 
band index b. The latter can be easily seen for s = 1 in the limit of infinite lattice. Replacing the sum in 


Eq. (19) by the integral and integrating by parts, we obtain 

. 1 dEbik) . „ 

Ji = -/ — r.—- sin(ka)dk . 

TT Jr, dk 


( 20 ) 


Since sin(A:a) > 0 in the integration interval, the sign of is determined by the derivative dEb{k)/dk which 
is positive for even b and negative for odd b. In the lowest Bloch band, is negative. Since this quantity 

plays an important role, we give it a special notation J = —JJi. 

Typical behavior of is shown in Fig. It is a decreasing function of Vb and s. Asymptotics of 
at large distances s is similar to that of the Wannier functions and given by 


17“ ~ |as| ^/^exp(-/io|as|) , 


( 21 ) 


where the constant Hq is the same as in Eq. (12). 

The dependence of J on the lattice amplitude Vb is shown in Fig. [^b). In the limit Vb —>■ 0, J can be 
calculated analytically using Eq. (11) and the result is 


lim - 

Vb-s-o Eyi 


( 22 ) 
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Figure 6: (color online) Tunneling matrix element in the lowest Bloch band, (a) The distance dependence for (from top to the 
bottom) Vo/F'r = 3 (red), 5 (green), 10 (blue), (b) Tunneling matrix element for the nearest lattice sites. Solid line is exact 
numerical result. Calculations within the Gaussian approximation [Eqs. ( |24| l] are shown by dashed-dotted line. Dashed line is 
the asymptotic formula ( |23| l. 


In the tight-binding limit Vb ^ which leads to Eq. (541, the width of the lowest Bloch band Eo(7r/a) — 
Eo{0) = 4J. Then from Eq. ^ we obtain asymptotic expression 


J 


( 9 



(23) 


As one can see in Fig. ib), Eq. ( |2^ gives quite accurate results for Vq/E^ > 10. 

In the Gaussian approximation (15) which is supposed to be valid for large Vq, the tunneling matrix 
element is given by 


J 

Er 


- 1 



(24) 


Although this equation describes correctly qualitative behavior, the magnitude of J appears to be underes¬ 
timated compared to the exact numerical results. 

Exact numerical results for J(Vo) were also fitted by the function with three parameters 


J 


= Pi 


Eft 


exp -p 3 



(25) 


which is motivated by Eqs. (23) and (24). For instance, in Ref. |123| we find pi = 1.39666, p 2 = 1.051, 
P 3 = 2.12104, and Ref. |124| suggests pi = 1.43, P 2 = 0.98, p^ = 2.07. Although both sets of parameters 
give indeed quite a high accuracy for Vb /Er > 3, the function (251 fails to reproduce the correct behavior 
for small Vq. Here we suggest a fit which has more parameters but works very well for small and for large 
Ho: 


J 

Er 


= Pi 


Efi 


exp 


-P3 


Eb. 


^exp 


-Pb 


V 

Eb) 


(26) 


with Pi = 0.116828, p 2 = 1.16938, p^ = 1.11717, p 4 = 0.63, ps = 0.369658, pe = 1.01448. 


2.2. Multi-dimensional lattices 


The generalization of the theory of one-dimensional lattices presented in Section 2.1 to arbitrary dimen¬ 
sion d can be done employing a standard approach from the solid state physics m\- Let xi denote the 
global minima of the periodic potential iy(x). We assume that the vectors Xi form a Bravais lattice and. 
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therefore, have the form xi = ^v^v-, where a.i, are primitive vectors that are in general not orthogonal 

to each other. The potential T^l(x) can be represented in the form of a Fourier series 


^l(x) = ^yjexp(igj -x) , 
j 


(27) 


where the coefficients Vj are given by 


= - / dxexp (-igj • x) Vl(x) . 


(28) 


gj’s in Eqs. (27 1, (28 1 are vectors of the reciprocal lattice determined by the conditions exp (igj • xi) = 1. In 
terms of the primitive vectors of the reciprocal lattice h^, defined by the identities 1 ^ 1^21 Sj 

have the following representation: gj = The integration in Eq. (281 is over one primitive cell of 

the volume v. 

According to the Bloch theorem, the wavefunction of the stationary Schrodinger equation has the form 
analogous to 


V'b(x; k) = U6(x; k)e' 


ik-: 


(x; k) = ^ ^ Cbj (k) exp (igj • x) . 


Ub 


(29) 


Imposing Born-von Karman boundary conditions + k) = k), we deduce that the wavevector 

k takes the values 

d 

kq = ^ , g. G z . (30) 

v—1 

The coefficients C{,j(k) and the energy eigenvalues Eb(k) are obtained from the solution of the eigenvalue 
problem 

^ (gj + k)^ C6j(k) + ^ 1/j/Cbj_j/(k) = Ef,(k)cbj(k) . (31) 

j' 


2M 


The Wannier functions can be determined in analogy to Eq. (§ as 


-1 


IF{,(x-xi) = I L,, j ^ •!/'&(x;k)( 


, —ik-xi 


(32) 




kelBZ 


with the orthonormality condition similar to Eq. (lOl. They allow to define the tunneling matrix 

In what follows, we restrict ourselves to hypercubic lattices that are created by d pairs of laser beams 
propagating along the axes. In order to create square or cubic lattices one has to avoid the interference 
of laser beams propagating in the orthogonal directions which is achieved if their frequencies are sufficiently 
different. This setup allows to create multi-dimensional lattices described by the potential 




(33) 


where d = 1,2,3 is the lattice dimension. The variables in the Schrodinger equation can be separated 
and the solutions are obtained from the one-dimensional ones according to the rules 


E = £'b(k) = '^Eb^iK) , ■0b(x;k) = tpb, {x^^K) 


(34) 
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atomic cloud 


Figure 7: (color online) Lin-6-lin laser configuration. 


Then the Wannier functions will be also given by the products 

d 

^bl(x) = (X^) . 


(35) 


Multidimensional analogues of Eqs. (18), (19) together with Eqs. (34), (35) allow to express the tunneling 
matrix in terms of the one-dimensional quantities as 


d 

ly—l 

which formally shows that do not vanish only along the lattice axes. 


(36) 


2.3. State-dependent potentials 

We consider two counterpropagating linearly polarized laser waves of equal amplitudes and frequencies 
with the wave number kj^, and the angle 0 (with 0 < 9 < 1 ^ 12 ) between the polarization vectors 

9 9 

e± = cos 2 ®2 ± sm - 63 . (37) 

The running laser waves form left- and right- polarized standing waves. This setup shown in Fig. is called 
lin- 6 *-lin laser configuration [55]. Here we consider a one-dimensional lattice but generalizations to higher 
dimensions are also possible |98L I126j . 

If the laser is tuned between the P 1/2 and P 3/2 electronic levels, the effective potentials acting on the 
ground-state sublevels of S 1/2 are different pil I127fll3()] . For instance. 


V\F= 2 ,mp=± 2 ){x) = V±{x) , (38) 

v=i.-.=±i>(^) = [^y±ix) + v^{x)]/4, 

where V±{x) = Vbcos^(fcLa; ± dl‘2). These potentials have the same period a = -njki^. For 0 = 0, they 
coincide but for other values of 9 their amplitudes are different and the positions of minima are shifted. This 
can be easily seen, if we rewrite V|i,±i)(a;) in the form 


V\F=i,m!,=±i) {x) = Vq cos^(fcLa; ± (f)) + Aq , 

^ Vl + 3cos^ 9 , (() = arctan tan 0 ^ , 


(39) 
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Figure 8: Scheme of the electronic transitions for lin-0-lin laser configuration. 


The band structure and all other single-particle states remain the same as in the case of large detuning 
discussed in section \2A\ 


2.4- Atoms with coupled ground states 

We consider again lin-0-lin laser configuration and turn to the case when the laser detuning A is com¬ 
parable to the hyperfine splitting of the electronic levels but still larger than the spontaneous emission rate. 
Assuming that only the ground states with F = 1 are populated and the laser frequency is close to the 
F = 1 ^ F' = 1 transition frequency of the D 1-line, the left- and right- polarized standing laser waves 
will couple internal ground and excited states with magnetic quantum numbers mp — 0, ±1 by V and A 
transitions, see Fig. 

The V and A laser-induced transitions lead to two sets of orthogonal Bloch eigenmodes which we denote 
by the indices 0 and A, respectively. The solutions of the Schrodinger equation are the three-component 
spinors of the form 

= (0, ^0,0)^, = (V'+, 0,, (40) 


which allow also to determine the Wannier spinors as well as the tunneling matrices for both types of modes 
according to Eqs. ([^, (181, (19l. 


2.4-1. 0-modes 

The effective potential acting on the atoms in the internal state a = 0 is given by 


Vb{x) = ^ L 


1 -|- cos 9 cos 


(-f) 


= tt/^l . 


(41) 


Therefore, these modes are exactly the same as those discussed in Sections |2.1.1[ |2.1.2[ provided that the 
amplitude Vq is replaced by Vq cos 9. 


2 . 4 . 2 . A-modes 

The non-vanishing components of the spinor ’ 5 'A) are solutions of the Schrodinger equation with the 
effective potential which has the form of a 2 x 2 matrix: 


VLix) = 


Vn 


nl 




f2±(a;) = cos ( tt— ± ^ 


(42) 


Due to the periodicity of the potential (421, they can be written down in the form of Eqs. with the 

scalars Ub and Cbn replaced by two-component vectors and , where the coefficients Cbn are solutions 
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-1 0 1-1 0 1-1 0 1-1 0 1-1 0 1 

fca/ir kaj'K kajir ka/n ka/ir 

Figure 9: Band structure of the A-modes. Vb = —8 -Er, 6 = 0° (a), 30° (b), 45° (c), 60° (d), 90° (e). 




fca/TT ka/iT kaj'K ka/K kajn 

Figure 10: Band structure of the A-modes. Vq = 800i?R, 6 = 0° (a), 30° (b), 45° (c), 60° (d), 90° (e). 


of the eigenvalue problem 


/ . ~ '-bn > 


(43) 


n ——oo 


= 

nn' 


p \ 1^0 / 1 cos 6» 

ifR(- + 2nj +-f ^ 


K) 


„i9 


1 ^ 1 1 A 

1 4" I 2 ^ giS ) On',n+l 


The lowest Bloch bands obtained by the numerical solution of Eq. (431 are shown in Figs. [9 10 (see also 
Refs. |131H133] b In contrast to the spinless case considered in section 2.1.1 the bands can overlap. There is 
a strong dependence on the angle 9 and the results are drastically different for positive and negative Vq. In 
the case of positive Vb the band gaps remain of the order of Er or vanish even for very large values of Vb in 
contrast to the case of negative Vb. 

The principal difference between the cases of positive and negative Vb can be understood if we apply the 
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Figure 11: Boundary between the regions of the normal and anomalous dispersion in the lowest Bloch band for the A-modes 
obtained from the condition = EQ^\n/a). In the shaded region, (Tr/a). 


unitary transformation 


= 


n 


-o_ 


fi = + 


(44) 


to the spinor ( 1 / 1 +, After the transformation we end up with the bright and dark states |182| . which 

are not degenerate in contrast to the original ones. The important point is that only the bright state is 
directly coupled to the electromagnetic field and influenced by the potential (411. 


We consider first the case ff — 0, when the transformation Z/f does not depend on the position z. In this 
case the Hamiltonian matrix is diagonal in the basis of bright and dark states and the dark state does 
not “feel" any periodic potential. Since for Vb > 0 the dark state has lower energy than the bright state the 
ground state is the same as for free atoms. If Vb < 0, the situation is reversed: The energy of the bright 
state is lower than that of the dark one and only the bright state is populated by the atoms. Therefore, 
increasing |Vb| one can strongly influence the lowest energy bands as in the case of spinless atoms. 

In the case 0 7 ^ 0, Z/f is a position-dependent transformation. The atomic center-of-mass motion leads to 
the gauge potential 

sin 6 > 


V, = En 


1 2 


(45) 


_ 1 -I- cos 0 cos (27rx/a) _ 

acting on the bright and dark atomic states and to the motional coupling of the states inn. The trans¬ 
formation Zi does not allow to diagonalize the Hamiltonian in the case 0^0. Nevertheless, it helps to 
understand what is going on, assuming that |Vb| S> i^R. In this approximation Vg <C |Vb|, and one can 
neglect the gauge potential for the bright state as well as the motional coupling between the bright and dark 


states m- Then the only potentials acting on the bright and dark states are given by Eqs. ( |4l| ) and ( [4^ , 
respectively. On the basis of the same argument as in the case 0 — 0 we see that the low-energy eigenstates 
in the cases Vb < 0 and Vb > 0 are determined by the potentials Vb and Vg, respectively. Accordingly in the 
case Vb < 0 the quantity Vb cos 0 defines the strength of the periodic potential, while in the opposite case 
the potential does not depend on Vq. 

This simplified description provides a correct physical insight but does not describe such an important 
feature of the A-modes as the change of the type of the dispersion relation in the lowest Bloch band 

under variation of the angle 0. In the case of spinless atoms, one always has a normal dispersion in the 
lowest Bloch band, i.e., dEo/dk > 0 for 0 < fc < tt/o. In the case we are dealing with, one can get anomalous 

H-e). The change of the dispersion 


dispersion, i.e., dE^^'^/dk < 0 for 0 < Zc < tt/o, as well, see Figs.l^e), 


type happens at the points (Vb, 0) indicated in Fig. 11 by the solid lines separating white and shaded areas. 
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The change of the dispersion types leads to the fact that the tunneling matrix elements for the nearest- 
neighboring sites in the lowest Bloch band 


Ja = - 




(46) 
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Figure 12: (color online) Tunneling matrix element for the A-modes in the lowest Bloch band for the nearest neighbors defined 
by Eq. l|46|l in the case of negative (a) and positive (b) Vo- 


can take positive and negative values. This can be seen from Eq. (19), which remains valid in the present 
situation as well, and it is demonstrated in Fig. by exact numerical calculations. For negative Vq and 
small 0, J\ Ri Jq but in general they are completely different. 


3. Spinless bosons with interactions 


We consider spin-polarized interacting bosons of the mass M in a periodic potential Vl(x). In all 
experiments with cold atoms in optical lattices there is also a (harmonic) trapping potential which can be 
of two different types. One of them is produced by the magneto-optical trap in three dimensions and has 
trapping frequencies of the order of 10 ... 100 Hz. This kind of trapping leads to the inhomogeneity of density 
profile of the atomic cloud which extends over several tens of the lattice periods. Another type of trapping 
can be produced by an additional optical lattice of large amplitude Vo_l in one or two dimensions which 
suppresses tunnelings in the corresponding directions and leads to effective trapping frequencies u}± up to 


~ 1 MHz as described by Eq. (141. The latter provides an opportunity to create (quasi) two-dimensional [8]- 
[T^ and (quasi) one-dimensional lattice systems and to reach experimentally the Tonks-Girardeau 

regime mm- 

In order to take into account all these possibilities, we will denote by Vt(x) the trapping potential of 
the first type and assume that the system’s dimension d can be less than three. Therefore, x is a vector in 
a d-dimensional space with d = 1, 2, 3. In the second quantization, the Hamiltonian has the form 


H = 


4'1'(x) 


-|^V2 + Vl(x) + Ft(x) 


'k(x) dx 


(47) 


J j ^^(x)4'^(x')Vat(x-x')^(x')^(x)dxdx', 


where 'I'(x) is a field operator which annihilates one atom at the point x. 
short-range, they can be described by two-body contact potential 

14t(x-x') =5d<5(x-x') , gd=- -, 53 = 

(a^V27r)' 


If the atomic interactions are 


M 


(48) 


with Os being s-wave scattering length in three dimensions and a_L the harmonic oscillator length corre¬ 
sponding to the frequency oj±. Eq. ( |48[ ) implies that the transverse confinement is not too tight, otherwise, 
confinement-induced resonances |135l 1135] may lead to modifications of the interaction parameter |137f - 
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3.1. Derivation of the Bose-Hubbard model 

The matter-field operator can be written down in terms of the Wannier functions as 


^(x) = ^ VTbi(x)abi , 


( 49 ) 


b,l 


where the annihilation and creation operators for the band b at site 1, Obi and aj^j, obey the bosonic 
commutation relations 


“bill 1 ojt, 1, 


= (5bib2^iii2 ) [abiii,ab2i2] = 0 ■ 


(50) 


Substituting (491 into the second-quantized Hamiltonian of interacting bosons in continuum (471, we obtain 


its discrete representation 

b lib 


E 


b2 
b 


+ 5 E E c'SIau"'''!,.,."! 


2b“b3l3“b4l4 ’ 


bib 2 lib 


bib2b3b4 libbb 


c/, 


bib2b3b4 _ 

llbbb “ 


H;b^= / W^biii(x)Hr(x)lTb2b(x')rfx, 

JJ fHb*iii(x)VFb* 2 b(^')Kt(x-x')hFb 3 i 3 (x')W^b 4 i 4 (x)dxdx' . 


(51) 


(52) 

(53) 


By doing different approximations in the Hamiltonian (511 one can derive lattice models of different types, 


such as multiband |140L I141| or extended |142| Bose-Hubbard models. These non-standard quantum lattice 
models were recently reviewed in Ref. |93| . 

The Hamiltonian can be simplified in the tight-binding limit Vq when the width of the Bloch 

bands becomes small and the gaps grow. In this regime we can keep only the terms corresponding to the 
lowest Bloch band b = 0. This is valid, if all the energy scales are less than the energy gap separating the 
first two Bloch bands. Since we will be dealing further only with the lowest Bloch band, we will drop the 
band index b. 

In this regime, the tunneling matrix elements with s > 2 become much smaller than fTf (see Fig. 
and, therefore, can be neglected. If the interactions are short-range, we can keep only local terms in the 


second part of Eq. (531. Since the spatial scale of the potential Vt (x) is much larger than the lattice period a 
and the Wannier functions are strongly localized, one can take Vt(x) in Eq. (521 out of the integral and use 
the orthonormality condition ( [Io| . In the isotropic cubic lattice, these approximations lead to the standard 
Bose-Hubbard Hamiltonian |47| 


Hbh = --^EE +h.c.) -f +E^T(xi)aj'a, , 


(54) 


v=l 1 


where e^, is a unit vector on the lattice in the direction v, J = is the tunneling matrix element for 

the nearest neighbors, Ud = is the on-site atom-atom interaction energy. For contact interaction (481 

it takes the form 


Ud = gd / |W^i(x)| dyi . 


(55) 


The tunneling parameter J was already discussed in section |2.1.3| The dependence of U on the lattice 
amplitude Vq is shown in Fig. 13 While J rapidly decreases, U grows. In the case Vq —>■ 0, Eq. (Ill yields 
analytical result for Ud- 


Ud f2Y 


3-d 

i-d f hu!± \ ^ Us 

^ ' ['Wj 7 


(56) 


Eqs. (22) and (561 show that the maximal value of the ratio J/Ud can be of the order of 10 in realistic 
experiments. Since J rapidly decreases with Vq and Ud grows, one can easily reach very small ratios of J /Ud 
which allows to access different regimes of the Bose-Hubbard model. 
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Figure 13: (a) On-site interaction constant in a three-dimensional lattice. Solid line is exact numerical result obtained according 

to Eq. ( |53[ |. Calculations within the Gaussian approximation [Eq. ( |57[ l] are shown by dashed-dotted line, (b) The ratio of the 
tunneling matrix element for the nearest neighbors in the lowest Bloch band to the on-site interaction constant. Solid line - 
exact numerical result, dashed line — estimations from Eqs. ( |23[ l, ( |57[ |, dashed-dotted line — Gaussian approximation | |24[ |, | |57[ l. 


In the Gaussian approximation (15) the interaction parameter takes the form 


Er 




Er ) \Er) a 


(57) 


The comparison with exact numerical results in Fig. shows that this expression overestimates the value 
of Ud, although it predicts correct qualitative behavior. 

Exact numerical results for the interaction parameter in three dimensions were fitted by 


Er \Er) a 


(58) 


with Pi = 2.985, p 2 = 0.88 |124| . However, this equation as well as Eq. (571 would imply that Ud vanishes 
in the limit Vq 0 which contradicts to Eq. (|56|). In order to find a fit which works well also for small Vq, 


we observe from Eqs. (551 and (561 that UdiVo) can be represented in the form 


cy 

Er 


8 

TT 


TT huj± 

4 



v^rA « 


(59) 


The numerical data are very accurately reproduced by a polynomial function u{x) = with the 

coefficients po = 8/27, pi = 0.554092, p 2 = 8.01432 x lO'^, P 3 = -8.94513 x lO'^, pi = 4.55577 x lO'^, 
P 5 = -1.12896 X 10-5, pe = 1.09512 x 10"^. 


3.2. Particle-number conservation 

Since the Bose-Hubbard Hamiltonian commutes with the operator of the total number of particles 


N = ^ , (60) 

i 

the latter is a good quantum number. Therefore, the eigenstates of the Hamiltonian can be represented as 

if the 

( 4 ) 


superpositions of the eigenstates of the operator (60 \ which are given by the product of Fock states 


|nrl = i 


yrrr/I 


■|0), T = 1,...,V, V = 


(7V + L- 1)! 
N\{L-1)\ ’ 


(61) 
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where T labels the configuration of the bosons and the occupation numbers of individual lattice sites nre 
satisfy the condition 

J2nri = N (62) 

i 


for any F. 


3.3. Translational invariance 

In the case of homogeneous lattices with periodic boundary conditions, the Hamiltonian commutes 
with the translation operator T which shifts the indices of the bosonic operators by one: •••—>■ 

• • •, and it holds ai = In order to simplify the notations, we consider here one¬ 

dimensional lattices but the formalism can be easily generalized to higher dimensions. Obviously, after L 
translations we always return to the original site, i.e., = I. Therefore, the eigenvalues of T are given 

by tk = exp(—liFa), where hK are the eigenvalues of the total momentum which takes discrete values 
K = Kq = determined by an integer q. 

The operator T commutes with the total-number operator N. The common eigenstates of T and N can 
be obtained acting by the projection operator [551 1143] 




j=0 


Tk 


on the states (611 which creates linear combinations of the form 

i/'r — 1 


|n_R-r) 


= H \.r 






j=0 


(63) 


(64) 


where t'r is a minimal number of translations required to map the state |nr) into itself which has to be a 
divider of L. Eq. (641 contains only those states |nr) that cannot be obtained from the others by cyclic 
permutations: |nr) ^ T^^ jnr') for j = 1,..., L — 1. The states (641 satisfy the orthonormality condition 
(nxr|nic'r') = Srr'^KK'- The sum over m in Eq. (64l implies that the states |ni<-r) 


do not necessarily exist 

for all values of Kq. For instance, the states with equal occupation numbers of all sites have j/p = 1 and, 
therefore, exist only for q = mL, m = 0, ±1,... On the other hand, the states with distinct occupation 


numbers have vy = L and, therefore, exist for any q. Eq. (641 can be also rewritten in the form 


nicr) = 



(65) 


Thus, the eigenstates of the homogeneous system under periodic boundary conditions can be labeled by 
two indices: K and H, where the second index H distinguishes between the states with the same value of K. 


The states (641 are used as a basis in the exact diagonalization because they allow to reduce the dimension 


of the Hamiltonian matrix by a factor of the order of the number of lattice sites. 


3 . 4 .. Momentum operators 

The translation operator T is a unitary operator and, therefore, can be represented in the form 

r = exp , (66) 

where is a component of the momentum operator which is a generator of translations on a discrete 
lattice in the direction v. Operators 11,^ should commute with H as well as N, and have the eigenvalues HK 
discussed in the previous section which are restricted modulo 2 'Kh/a. 
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We consider first the usual momentum operator in a continuous space 


P = y ^^(x) (—iSV)'I'(x) dx . 


( 67 ) 


A lattice version of P can be obtained using the basis of Wannier functions. In the tight-binding approxi¬ 
mation, we obtain 


where 


P = iPo ^ e,, ^ («l'+e,.«i - h-c.) , 
1 


Pq = h / W{x) — W{x — a) dx . 

ux 


In the Fourier space introduced via the transformation 

1 


27r 


ak = ^‘/?ki“i> <Pki = exp (ik • xi) , > 


the momentum operator P takes the following form 


( 68 ) 


(69) 


P = 2Po OkOk sin {Ka) ■ 


(70) 


U = 1 


Note that the operators in the momentum space hk satisfy standard bosonic commutation relations and 
have the following property 

— 




kelBZ 


Operator P cannot be identified with II, because it does not commute with the Hamiltonian (54 \ due 
to the interaction term. Instead, the quasi-momentum operator 


Q = fikd 


k“k ; 


kGlBZ 


is introduced which commutes with the Hamiltonian as well as with the translation operator. However, it 
is not restricted modulo 2TTh/a. A proper momentum operator which satisfies all the conditions is given 

by [ml El] 




L-1 


• IJ 


exp - 1 


(72) 


Although being of fundamental importance, the explicit form of the momentum operator does not play a 
role for the interpretation of the experiments with ultracold atoms in optical lattices. What is more relevant 
is the quasi-momentum distribution |145| determined as 


1 


P(k) =-(^t(k)^(k)) ^ ^(k) = 


1 


'k(x) exp (—*k • x) dx 


(73) 


where N should be understood in general as an expectation value of the operator N. Taking into account 


Eqs. (491 and (691 it can be written down in the form 


P(k)= W{k) 


^ (Qlc«k) 

N 


^P(k) = l 


(74) 
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where 


W (k) = —^ / W (x) exp (—ik ■ x) dx , 

a^i z I 


( 75 ) 


and the values of k are not restricted to the first Brillouin zone. The normalization constants in Eqs. 
are chosen such that 


^4^(k)§(k)= / ^^(x)^(x)dx = iV . 


In the Gaussian approximation (151, Wik) has the form 

a 


W{k) = (Anf* 


exp ( --k ttho 


(76) 


Usually in the theoretical analysis of the lattice problems based on the Bose-Hubbard model (54), the 


Wannier functions are omitted. Then the expression for the quasi-momentum distribution is considered to 
be 


P(k) = 


N 
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p(k) = 1. 


(77) 


In an infinite system, k is a continuous variable and it makes sense to redefine the quasi-momentum distri¬ 
bution as Poo(k) = (Pa/27r)‘^P(k) in order to have a proper normalization: -Poo(k)(ik = 1. 

Eq. (691 allows to express the momentum distribution in terms of the operators in real space in the form 


of a double sum over the lattice sites. In a translationally invariant system, the double sum can be converted 
into a single sum and Eq. takes the form of the discrete Fourier transform: 


L-l 


L-l 


= JV ^ 

fi=0 


exp [zkq • (xi - xo)] 


^0^1 / 


(78) 


id=0 


One can show that P(kq) in the last equation is real-valued. 

In general, the quasi-momentum distribution Poo (k) defined by Eq. (771 is an even and periodic function 
of ku with the period 27: ja. For J > 0, it takes maximal (minimal) values at ka = ttiti, m = 0, ±2, ±4,... 
(m = ±1, ±3,...). In the case of neg ativ e J, the positions of the minima and maxima are reversed. How¬ 
ever, the presence of |IU(k)p in Eq. (74) for the true quasi-momentum distribution destroys this periodic 


structure resulting in smaller heights of the peaks with larger values of k. With the decrease of the hop¬ 
ping parameter J the spatial correlations of bosons become weaker which leads to the broadening of the 
momentum distribution. 


4. Basic definitions 

4 ..I. Thermodynamic quantities 

We remind the definitions of the basic thermodynamic quantities which are often used in the studies of 
the many-body lattice problems. Starting with the grand-canonical partition function 


Z(^,T)=Trexp - 


H- nN 
knT 


the free energy is defined as 

Pifi,T) = -kBTlnZifi,T) . 

The derivative of the free energy gives the mean number of particles per lattice site (filling factor) 


dP 


N 


(hi) = -L — , (hi) = , 


dyi 
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(79) 

(80) 

(81) 









and the derivative of the latter gives the (isothermal) compressibility 

■ 

\ * 9 /^ / T 

The compressibility is related to the fluctuations of the total number of particles: 

1 (iV2) _ iV2 

~ 1 ? ^ ■ 

At zero temperature, the free energy takes the form 

= 0) = min (Ejsf — fiN) , 

N 


( 82 ) 


(83) 


(84) 


where Aat is the ground-state energy of N particles with N being a non-negative integer. If N can be 

(lization : 

/r(iV) = 


considered as a continuous variable, the minimization in Eq. (841 gives 

dEjq 


dN ■ 


(85) 


However, strictly speaking N in Eq. (|84| is discrete and the minimization leads to the fact that the chemical 
potential /r is enclosed in the interval [n-, /r+], where the boundaries are given by 


± {En±i — Em) ■ 


Discrete analogue of Eq. (821 has the form 

K-i = [/r+(fV) - 

and defines the one-particle (’charge’) gap 

Ac = Em+1 + Em-1 — ‘2Em ■ 


( 86 ) 


(87) 


( 88 ) 


This is different from the ‘neutral’ gap A„ which is defined as the energy difference of the lowest excited 
state and the ground state with the same number of particles N. 

In inhomogeneous lattices, it is useful to consider local quantities. For instance, local (on-site) compress¬ 
ibility [TTzlITIH] 

(89) 


Ki = 




quantifies the response of local density to the local variations of the chemical potential. 


4-2. Superfluidity 

Superfluidity is one of the most fascinating phenomena which can occur in many-body quantum systems 
at low temperature |1491 1150] . It was discovered first in the experiments with liquid helium |151L 1152] 
and also observed with ultracold atomic gases in traps |153H156] and in shallow optical lattices |157fll59] . 
Superfluidity is usually related to the flow properties of a quantum system which is assumed to be composed 
of a normal and superfluid components distinguished through their behavior in the presence of moving 
boundaries. If the system is enclosed in a narrow region between two moving walls, the normal component 


^ This is slightly different from the standard definition [US] 

__1 fdV\ _ V fdN\ _ a<^ fd{ni)\ _ a<^ 

V \dpJj,jy ydpLjj^y (ni}2 V dp. ) j, y (rii}2 
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is dragged by the walls, whereas the superfluid remains at rest. In other words, the superfluid is at rest in 
the lab frame, while the normal component is at rest in the frame of the moving walls. 

The superfluidity is quantified imposing twisted boundary conditions on the many-body wavefunc- 
tion 

4'(... ,Xj + L'e^, ■■■) = .. ,x^-,...) , (90) 

where L' is the linear size of the systeirj^and the twist angle 9 € (0, tt). This requirement leads to the increase 
of the free energy T which is attributed to the kinetic energy of the superfluid. Since the corresponding 
velocity is fixed by the value of theta, the number of particles in the superfluid component is determined 
as 


J-(k,) - J-(O) = 


■-N, 


ks — 


(91) 


2M ’ L' ’ 

which readily gives the superfluid fraction determined as i/g = N^/N. According to this definition, Ug G [0,1] 
and the upper bound corresponds to noninteracting particles in free space at zero temperature. In the limit 
kg —t 0, Eq. (91 \ yields 


Then in the case of noninteracting particles in a periodic potential at T = 0, we get Ug = M/M* |16311164] . 
where M* is the effective mass, which follows from the dispersion relation of non-interacting particles. In 
the tight-binding regime, M* is given by 


M* = , 


(93) 


and Vg = TT^ J/Er is exponentially small. This makes the measurement of the superfluid fraction in deep 
optical lattices rather difficult. Nevertheless, this quantity plays a fundamental role in the theoretical studies 
due to the possibility of the superfluid-insulator transition. 

In a lattice model, the dispersion relation of noninteracting particles Ck has a quadratic form only for 
small k. For the description of the superfluid properties of the system it is convenient to introduce the 
superfluid stiffness in analogy to Eq. 


fh = 


J^jkg) - J-(O) 
N (ck, - eo) 


k -e A 
Kg — ^1/ r 

La 


(94) 


For noninteracting particles in a translationally invariant lattice, = 1 for arbitrary kg. For small kg 


Eq. (941 can be rewritten in the form (921 with M replaced by the effective mass M*. In the hydrodynamic 


approach the superfluid stiffness is related to isothermal compressibility k as 

f% = M^Kcl/{ni) , 


(95) 


where Cg is the sound velocity. 

From the practical point of view, twisted boundary conditions can be introduced in slightly different ways. 
One of the often used possibilities is to impose antiperiodic boundary conditions |165| which corresponds 
to 0 = TT. This allows to avoid complex arithmetics in the numerical calculations. Another option is 


to introduce Peierls phase factors exp(i6*/L) into the hopping part of the Hamiltonian (541 by means of 
transformation di —>■ a\ exp (ikg • xi). It is interesting to note that the Peierls phases of arbitrary magnitude 
can be created by periodic shaking of the optical lattice along a chosen direction m as it was demonstrated 
experimentally in Ref. nszl. 

One can also consider the limit 0 —>■ 0 which leads to the concept of winding number |168] used in 
QMC calculations. In addition, in this case one can express the superfluid stiffness at zero temperature 
as |B51 [Ml IT55] 

|(^A 


fN = - 


2NJ 




H. 


kin 


A^O 


J. 


4'n 


E\ — En 


(96) 


^In the case of a periodic potential, L' = La. 
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where 


Ji, = iJ 


1 


(97) 


is the current operator for the direction v which up to a constant prefactor coincides with the operator ( 681 . 


The operator can be formally obtained taking the derivative with respect to 6 in the Hamiltonian with 


the Peierls phase factors and considering the limit 0 —0 (see, e.g., Ref. |170| b Eq. (96l shows that the 
superfluidity is not just a property of the ground state but also contains information about all excited states 
described by the second term which always gives a negative contribution. 

4-.3. Bose-Einstein condensation 

The Bose-Einstein condensation is in general defined as a macroscopic population of an eigenstate of the 
one-body density matrix |171LI172] /9i(x,x') = ('iE'^(x)^(x')). Formally, this is determined from the solution 
of the eigenvalue problem 

y Pi(x,x')<?i'*(x')dx' = , (98) 

where index i labels the eigenstates. 

In the tight-binding regime, the eigenfunctions (()i(x) as well as the held operator 'I'(x) are decomposed 


in the basis of the Wannier functions for the lowest Bloch band which leads to the lattice version of Eq. (981: 


1 ' 


(99) 


where the eigenvalues Ni remain the same. The discrete eigenfunctions satisfy the orthonormality condition 


1 




( 100 ) 


and the sum of all Ni gives the total number of particles N. If the largest eigenvalue is labeled by i = 0, 
the condensate fraction is determined as = Nq/N. In a homogeneous lattice with periodic boundary 
conditions, the discrete one-body density matrix with the entries ^^(li, I 2 ) = depends only on I 2 — R. 

The eigenfunction corresponding to the largest eigenvalue is constant, (poi = and 


iVn = 


L-1 

E- 

^1=0 


L-l 

•E< 

id=o 


) 


(101) 


Therefore, /^ = P(0). 

A sufflcient condition for the existence of a Bose Einstein condensate is the off-diagonal long-range order 
of the one-body density matrix [891117311174) . If the asymptotic value 




0 = 1™ Pi(x,x') = A^P(k = 0) 

|x—x' |—>-oo 


( 102 ) 


does not vanish, there is a finite fraction of particles with zero momentum. In the tight-binding regime and 
in a translationally invariant lattice, 


N^=°/N= W{0) 


/: 


N 


(103) 


Since 


1T(0) 


is smaller than one, the fraction of particles with zero momentum is smaller than the con¬ 


densate fraction in agreement with a general statement that is a lower estimate of Nq ina- It was 

also demonstrated in two dimensional lattices without using tight-binding approximation that the difference 
between and can be large |175| . In the experiments with ultracold atoms, P(k) is usually 

measured which gives a direct access to (see, e.g., |22l[2!?l[2Hlll76p . How ever, due to inhomogeneities 

caused, for instance, by harmonic confinement, there is no simple relation like (1031 between and 


and more careful analysis is necessary in order to extract from the experimental data. 


26 






















5. Main experimental techniques 

Possibility of experimental control is an attractive feature of ultracold atoms in optical lattices. Below 
we briefly discuss the main experimental methods which provide information about the state of the system 
and give motivation for theoretical studies. 


5.1. Time-of-flight imaging 

Spatial coherence of ultracold atoms is usually probed by the time-of-flight imaging of expanding atomic 
cloud after sudden switching-off the lattice and confining potentials |51I^[T771. The images obtained by 
the resonant absorption of photons are directly related to the density profiles. When the phase coherence 
length is at least of the order of several lattice spacings, the density distribution of an expanding cloud shows 
an interference pattern which has the symmetry of the reciprocal lattice. This is usually interpreted as a 
signature of superfluidity. If the phase coherence length is of the order of one lattice spacing, the density 
distribution becomes broad and does not show any peaks 1371 il]. 

In the experiments, the expansion times are sufficiently long such that the interaction effects during the 
expansion are negligible |178| . For 


o-JRt > —\ - - , 

a V hiOiio 


(104) 


where Rq is the characteristic size of the atomic cloud before expansion, and Who is the effective oscillation 


frequency at the bottom of a lattice well given by Eq. (141, the matter-field operator takes the form |178II179] 




M\ 


d/2 


— j W(k);^exp(-fk.xi + f—x: 


M 


2hi 


^1 ’ 


k = 


Mx 

ht 


(105) 


where we omitted unimportant phase factors. The expectation value of the density operator (i/ji (x, t)'0(x, t)) = 
Ptof(x, t) is given by 


fM\‘^ 


u 

lF(k) 


I1.I2 


ik ■ (xp - Xb) - (xp - xf. 


ii«b) ■ 


The second term of the exponential function in Eq. (1071 can be neglected, if 


WRt > 


icRo TT^ 


(106) 

(107) 

(108) 


where £c is a characteristic coherence length of the system which is of the order of fe w a in the Mott-insulator 
(MI) phase and of the order of Rq in the superfluid (SF) phase. Under condition (1081, 5(k) = 
and, therefore, describes the quasi-momentum distribution However, the expansion times in the 

experiments are usually shorter than those given by the “far-fleld" condition (108), although Eq. (104) is 


fulfilled. Therefore, for the interpretation of the experimental data it is important to keep the second term 


in the exponential function in Eq. (1071 as it was demonstrated in Ref. It is also necessary to keep in 

mind that experimentally one observes a two-dimensional column density p±{'r±,t) obtained from pxoF(r, t) 
by the integration along the probe line of sight. 

The structure of the density distribution of the time-of-flight images is quantitatively described by 
visibility |178II18()| 

P±(kmax)-P±(kmin) _ 

P± (kniax ) +P-L(k min ) 

In a three-dimensional setup, a special choice of the two-dimensio nal v ectors k^ax and kmin allows to cancel 
the contribution of the function Wo(k) and replace p_L(k) in Eq. (109) by ^^(k). This is achieved, provided 


that kmax is iu the center of the second Brillouin zone, i.e., k^axO = (27r,0), and kj^jn is along a diagonal 
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and has the same length as kmax, he., kminQ = ■\/2(7r,7r). Eqs. (1061, (1071, (1091 show that the one-body 
density matrix plays an important role in the experiments with ultracold atoms in optical lattices. 

An interesting aspect of the time-of-flight images is that in each experimental measurement a single 
realization of the density distribution is observed rather than the expectation value |181| . This allows to 
extract the density-density correlation function 


0(xi,X2) = ('0'l'(xi,t)?^(xi,t)?^'l'(x2,t)-!/>(x2,t)) - (xi, t)?/;(xi, t)) ('01'(x2, t)'0(x2, t)) 


( 110 ) 


from the experimental data and constitutes the basic idea of noise-correlation interferometry |182fll84] . It 
provides information on spatial order in the lattice that is absent in the average density. 

In order to simplify equations, we consider the regime of large expansion times such that the condi¬ 
tion (1081 is fulfilled. Using equal-time commutation relations for bosonic field operators and Eq. (105), we 
can rewrite Eq. (1101 in the form |184l I185| 




2 [ 


u) 

lU(ki) 

I 

lU(k2) 


■ E 

I 1 I 2 I 3 I 4 


:‘1-X'3)e*k= (x‘2-X.4)(at al a, a, 

\ li I 2 I 3 I 4 


+ (5(ki-k2)^e*k4.(x.,-x.,)^.t^-^^^_ ^(k2) "5(ki)5(k2) 


1112 


( 111 ) 


where the first term is the second-order correlator, while the second term which gives delta-peak for vanishing 
relative momentum corresponds to the correlations of an atom with itself (autocorrelation). In the case of 
fermions, this term would enter with the minus sign. 


In terms of the operators (69) in the quasi-momentum space, we get 

2d 

W(ki 






( 112 ) 


X 




W{k2) - («ki)(nk2) - (^ki)<5k,-k2.q|^ + («ki)(5ki 


where q is a d-dimensional vector of arbitrary integers. This leads to the following lattice version of the 
noise-correlation function used in theoretical studies |185L 1186] : 

0L(ki,k2) = (hkihk^) - (nki)(hk2) - (nki)^ki-k2,q|;^ > q 7^ 0 . (113) 

If the lattice potential is switched off slowly enough such that no transitions between the Bloch bands 
take place, the quasi-momentum of atoms is projected to the usual momentum. As a consequence, the 
population of the nth Bloch band is mapped into the momentum interval corresponding to the nth Brillouin 
zone. This is used in the bandmapping technique which allows momentum-resolved measurements of the 
populations of the lowest Bloch bands [3711187| . 


5.2. Optical Bragg spectroscopy 

Information about the excitation spectrum of the many-body system can be obtained with the aid of 
optical Bragg spectroscopy fKlimRT UTTniTTTl II8811191] . In ultracold atomic gases it is performed using 
two laser beams with the wavevectors ki and k 2 and frequencies uji and UJ 2 . This setup allows to transfer 
the momentum Kk = h{ki — k 2 ) and the energy hui = h(u}i — 0 J 2 ) to the atomic sample which can be tuned 
independently |40| . The measurement is performed again as a time-of-flight imaging of the expanding atomic 
cloud. 

This kind of perturbation is described by the Hamiltonian 


.^Bragg = Ueragg J dx COS (k • X - Ujt) (x, t)tp{x, t) 

_ ^Bragg [ — iujt) + p{k, t) exp {iujt) 


( 114 ) 
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where 


p(k, t) = / dxexp (—ik • x)'0’^(x, t)^(x, t) , p^(k, t) = p(—k, t) . 


(115) 


In the linear-response regime, the fluctuation of the density induced by the perturbation is given by the 
susceptibility x(k, cu) which depends only on the properties of the system in the absence of perturbation. It 
is determined by the relation m 


(p(k,t))vB.agg - (p(k,t)}o = [x(k,w)e -h x(k, Vb 


ragg • 


(116) 


Due to the causality of the response to the perturbation, the susceptibility x(k, w) is an analytic function 
of uj in the upper half of the complex plane and, therefore, satisfies the Kramers-Kronig relation 


1 V 

X(k,a;) = — / x(k,a;')^- dw', 

in uj' -UJ 


(117) 


where V denotes the principal value. Eq. (117) establishes the relation between the real and imaginary parts 
of X- 

The probability to transfer the momentum Kk and energy huj into the many-body system is proportional 
to the dynamic structure factor 


s(k,-) = T 


dt (Ap(k, 0)Ap(—k, t)) exp(—iwt) 


where Ap(k, t) is the spatial Fourier transform of the density-fluctuation operator 

Ap{x,t) = - {tp^x,t)'>p{x,t)} . 

In terms of S(k,uj), the susceptibility function is given by 


(118) 


(119) 


Re[x(k,a;)] = f 

J — ( 


S{k,uj')^’^ + S{-k,uj')-^ 


UJ' — iU 

Im[x(k,a;)] = tt [S'(k, w) - 5'(-k,-w)] . 


uj' + UJ 


duj' 


( 120 ) 


If we denote the eigenstates of the unperturbed Hamiltonian by I'I'a) with A = 0 corresponding to the ground 

in the form 

d{uj - ujxo) ■ 


state, the expression (118) at zero temperature can be rewritten in the form 

5(k,c.) = J^|(vI/A|A^(-k)|vI/ 


^0/ 


A 


A finite-temperature generalization of Eq. (121) is given by 


5(k,w) = ^ ^exp(- 


hZ 


X1X2 


E\2 

knT 


(^Ai|Ap(-k)|TA2) (5 (w-waiA2) . 


( 121 ) 


( 122 ) 


where Z is the partition function. 

The dynamic structure factor obeys certain sum-rules |1711 1192] . For instance, the energy-weighted 
moment is given by 


^(k,^)^^^ = -(^Ap1'(k), Ap(k)jj) = , 


(123) 


which is known as f-sum rule. It holds for a wide class of many-body systems independent on the external 
potential and temperature. Eq. (123) is directly related to the equation of continuity and the conservation 
of the particles number. The inverse energy-weighted moment in the long-wave length limit has the form 


lim 

k^O 


UJ 2 


( 124 ) 
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where kt is the isothermal compressibility (see footnote in Sec. 4.1 1 . Eq. (1241 is referred to as the com¬ 
pressibility sum rule. The zeroth moment of S'(k, w) yields the static structure factor: 


^o(k) = I 


= ^(Ap(k)Ap(-k)) 


(125) 


which is the Fourier transform of the density-density correlation function in real space. Together with 
Eq. (1151, this implies 

S'o(O) = ((iV^) - /N . (126) 

If the fluctuations of the total number of particles vanish, which is always the case in the canonical ensemble, 
S'o(O) vanishes as well. 

If we restrict ourselves to the lowest Bloch band, the static structure factor takes the form 


5o(k) = l + ^ ^ Gq,,(k)Gr,i3(k) 


N 


{ai^ a,ya,^a,,) — (ai_ a,^)(aya, 


U-i U-i U/I 

h I 3 h I 4 




(127) 


I 1 I 2 I 3 I 4 


where 


G 1 P 2 (k) = J W exp (-tk • x) (128) 

with ITi(x) = IE(x — xi) being the Wannier function for the lowest Bloch band. Thus, S'o(k) contains 
correlations of four lattice points. It is easy to see that Gqi 2 (k = 0) = and Eq. (1261 is fulfilled. 

The absolute values of Gqi 2 (k) decrease with increasing distance between the lattice points li and I 2 . 
Taking into account only the dominant terms with I 2 = li and I 4 = I 3 in Eq. (1271, we come to the expression 


where 


5o(k) = 1 + Gg(k) ^o(k)-l 


-5o(k) = X! I 2 ) exp [ik • (xq - xq)] 


I1.I2 


with the particle-number correlation function 

Fn{h,h) = (nqhq) - (nq)(nq) 


(129) 


(130) 


(131) 


is the discrete analogue of S'o(k). It corresponds to the fluctuations of the particle number described by the 
operator Ahi = fii — (hi) and contains only particle-number correlations of two sites. 

Go(k) in Eq. (1291 is defined as 


Go(k) = |Gii(k)| = J dx|IE(x)|^exp( 


—zk • x) 


In the Gaussian approximation (151 it takes the form 

Go(k) « exp 


fca V 

^ ) 



(132) 


(133) 


The “discrete" d ynam ic stru c ture factor S'(k,w) corresponding to the static structure factor 5o(k) is 
determined by Eqs. (118 1 , (121 1 , (1221 with the operator Ap(k, t) replaced by 


Ah(k) = (hi — (hi)) exp (—zk ■ 


XI 


(134) 
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This corresponds to the perturbation described by the Hamiltonian 

-^Bragg — ^^ragg / ^ COS (k • Xj Ojt') 7i\ = 




(k)e-“* + h(k)e" 


(135) 


The sum rule (1251 is also fulfilled for S'(k, w) and iS'o(k): 


h 


^(k,a;)da; = — (An(k)An(—k)) , 


(136) 


and Eq. (126 \ holds for S'o(k = 0) too. However, the f-sum rule ( |123[ ) and the compressibility sum rule (124 \ 
take for S'(k,a;) slightly different forms. The f-sum rule becomes |65| 

r°° ~ 1 rr - .. n . n 

S{k,u})u}du} = -( An(k),HBH ,An^(k)) 




kinx 


i/=l 


where 


= - J 


E(4 




(137) 

(138) 


is the kinetic-energy part of the Hamiltonian for the direction v. In a homogeneous isotropic lattice, Eq. (1371 
can be rewritten in the form 


I S(k,uj)ujduj = 


where 




(139) 


(140) 


is the energy of free particles. Comparing with Eq. ( |123 1 we see that the single-particle dispersion relation 
in continuous space h?k‘^/{2M) is replaced by the corresponding tight-binding dispersion relation ek and 
instead of the total particle number N we have (ajai_|_e )L‘^ which is less than N. The compressibility sum 
rule in a lattice has the form 


lim 

k^-O 


LU 2 


(141) 


where k is defined by Eq. (82). 

In some rather general cases, the dynamic structure factor is determined by a single excitation mode 
with the energy Huj^. Then at T = 0 it can be approximated as 


^(k, w) = 2k<J(w — Wk) . 


Using the sum rules (|136|), (|139|) we obtain a lattice analogue of the Feynman relation 

Ck 




hojk = 


) 


(ni) ^o(k) 


Employing in addition the compressibility sum rule (141) we obtain 


lim Wk = Cs |k| 
k-s-O 


Cs = X 




n\i K 

where Cg is the sound velocity. From Eqs. (1431, ( 144[ ) we get 




(142) 


(143) 


(144) 


(145) 


For small |k|, G'o(k) can be decomposed in powers of |k| . Therefore, the behavior of S'o(k) in the limit 
k —)■ 0 is also described by Eq. (1451. 
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5.3. In-situ imaging 

If the absorption images are taken without time-of-flight, they reveal in situ atomic density distribution [51 
m]. This technique allows to achieve the spatial resolution of the order of few microns which makes possible 
to resolve individual sites only in some special cases [Ml I193| . Below we give a brief overview of other 
methods that reach higher resolution and in addition allow efficient manipulation on a single-atom level. 
These methods allow to access the particle-number statistics of the individual sites iHHiniiii] as well as 
nonlocal correlation functions |31[ 119411195] . 

5.3.1. Microwave spectroscopy 

Microwave spectroscopy is based on the resonant transfer of atoms from one internal state to the other 
with the aid of two-photon pulses composed of microwave and radiofrequency photons. In the experiment of 
Ref. |196| . these were the states |1) = |F = l,mF = —1) and |2) = |F = 2,mF = 1) of ®^Rb. In the presence 
of interactions, the transition frequency is shifted by Aw ~ p(x)(a 2 i — an), where p(x) is the density, 021 
is the scattering length of two atoms in the states |1) and |2), an is the scattering length of two atoms in 
the state |1). After the pulse, the atoms in the state |2) are detected by light absorption imaging showing 
the spatial distribution of the atoms with a given density. Changing the photon frequency and repeating 
the measurements one can reconstruct the complete density profile. 

5.3.2. Scanning electron microscopy 

The fundamental physical process of the scanning electron microscopy is ionization of atoms by a focused 
electron beam. The ionized atoms are extracted from the atomic sample with an electrostatic field and 
subsequently registered by an ion detector. The diameter of the electron beam determines the spatial 
resolution which was 100 — 150 nm in the experiments of Refs. fWl [Ml. The ion detection provides 
single atom sensitivity and allows to count atoms on individual sites. This method can be used not only for 
measurements of the atomic distribution but also to produce arbitrary patterns of occupied lattice sites |198| . 

5.3.3. Fluorescence imaging 

Fluorescence imaging is a technique that allows to measure occupation numbers of the lattice sites 
modulo two. The atomic sample is illuminated with a near-resonant light which provides simultaneously 
sub-Doppler cooling pmii^iTM]. During this process, pairs of atoms undergo light assisted collisions 
and quickly leave the trap before they can be detected. After that only single atoms remain on the lattice 
sites with odd initial populations and the sites with initially even populations are empty. The remaining 
atoms scatter several thousand photons during the exposure time and can be detected with high fidelity El- 
Spatial resolution of the obtained images is about 600 — 700 nm mmm- 

The physical observable measured by the fluorescence imaging is described by the parity operator 

Si = exp(i7rhi) = (—1)"' (146) 

which yields —1 and -|-1 for odd and even occupation numbers, respectively. Measuring the parities at 
different lattice sites and averaging over many experimental realizations allows to determine the parity 
correlation function |199| 

(li,l2) = (spsij) - (sp)(si2) . (147) 

It was measured in one- and two-dimensional lattices |1941 I195j and provides an important information 
about correlated particle-hole pairs. 

The limitations imposed by the parity projection were circumvented in bilayer Bose-Hubbard systems. 
Using the interaction blockade in double wells [MTI . the occupation-sensitive transport between the two 
lattice planes was engineered that made possible to resolve the occupation numbers of lattice sites in one 
plane n\ = 0,1, 2, 3 and to observe the spin ordering in a mixture of two hyperfine states |201| . 

6. Simple special cases 

In this section, we discuss limiting cases which allow exact analytical solutions or do not require much 
computational effort. 


32 






















6.1. Ideal Bose gas 

We start the discussion of the many-body phenomena considering first an ideal gas of N bosons in a 
finite lattice of L‘^ sites under periodic boundary conditions (ai+e„L = oi)- In this non-interacting limit 


([/ = 0), the system is governed by the first term of the Hamiltonian (54). In the homogeneous lattice, it is 
convenient to work in the momentum representation ( |69[ ), where the Hamiltonian takes the diagonal form 

(148) 

(149) 

(150) 


H = J2 . 

k 

which readily gives the energy eigenvalues of a single particle 

d 

Ck = —2 J cos (fcyu) 




and the corresponding eigenstates 


|k) = aj^|0) . 


6.1.1. Energy spectrum 

From the solution of the single-particle eigenvalue problem, one can construct the eigenstates of the 
many-body system as products of the Fock states 

|h) = (^|hk), yhk = iV, (151) 

k k 

with the occupation numbers hk of the k-modes. They have the energies 

Eh = . (152) 

k 

In the ground state, all the bosons occupy the single-particle mode with k = 0. In the lowest excited state, 
N —1 bosons occupy the mode with k = 0 and one boson is in the mode with k^a = Isil'n jL in any direction 
V. Thus, the lowest excited state is 2d-fold degenerate and has the total momentum ±2Trh/(La). If L 
is large, all the eigenenergies are infinitesimally close to each other and, therefore, there is no gap in the 
excitation spectrum in the thermodynamic limit. Note that the energy of particle-hole excitations vanishes 
for any lattice size. In general, the energy spectrum has a band structure and its explicit form depends on 
the number of particles. In order to demonstrate this, we consider two examples. For the sake of simplicity, 
we shall restrict ourselves to one dimension but the generalization to higher dimensions is straightforward. 

In the case of two atoms, the energy spectrum is obtained by addition of contributions from two free 
particles with the momenta P± = ± p), where P = P+ + P_ = hK is the center-of-mass momentum and 

p = (P+ — P_)/2 = hk is the relative momentum. The energies are given by 

= EK,k=eK_^+€K_^i^ = -qK cos {ka) , qk = 4:Jcos , (153) 


and form a quasi-continuous band with the boundaries Ex^±-Kja = 

In the case of a commensurate filling, N = nP, with n being an integer, the lower boundary of the band 
is given by 


E. 


K,min _ 

nL 


2JN 

\ sinTn 


\L/ irn 


< Ka < 


(154) 


and the upper one 


TP K max. 
EjiL — 


—E^^™ , if P is even , 

2 JN cos (j) , if P is odd . 


(155) 


In the thermodynamic limit, all the eigenstates of the 7V-body system are enclosed in the interval [—2 JiV, 2JN] 
for any value of K. 
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6.1.2. Ground-state properties 

The ground state of N ideal bosons 


l^o) = |A^4.o) = 

k 




-| 0 ) 


(156) 


is unique and has the energy En^ = —2dJN. It is a perfect superfluid state with = 1. Taking into 

account Eq. (691, the expression (1561 for the ground state can be rewritten in terms of the occupation 
numbers of the individual lattice sites ni as 

|no> = ^ 


Nl 


\n\) 


V vW 

n 1 

Therefore, the number of bosons n\ at any lattice site 1 follows the binomial probability distribution 


(157) 


p[n\ = n) = 


m 


1 


i\(N — n)\ \L‘^ 


1 - 


1 

I? 


N-r. 


n = 0, l,...,iV, 


with the mean value (hi) = N/L'^ and standard deviation 


= y {n^i) - {niV = ]J (ni) ■ 


(158) 


(159) 


Eq. (1581 reflects the fact that in the absence of interactions every boson can be placed on a particular 
site independent of the others with the probability l/L'^. In the thermodynamic limit, the probability 


distribution (1581 takes the form 


p(ni = n) = e 


{n\y 


n\ 


This result is known as Poisson limit theorem. 
The one-body density matrix has the entries 




(160) 


(161) 


for any li and I 2 and does not show any dependence on the system size and dimensionality. Therefore, the 
condensate fraction is exactly one. 

The particle-number correlation function has the following form: 


Fn (ll,l2) = (hi) ( 5li,l2 - 


(162) 


For li = I 2 , Eq. (162) is consistent with (159). The term l/L'^ makes no contribution in the thermodynamic 
limit and F„(li,l 2 yf li) vanishes. However, it gives singular contributions to the static structure factor. 
From Eqs. (129) and (130) we obtain 


,So(k) = l-5k,^q, 5o(k) = l-G2(k)4,2, 




(163) 


,So(k) and 5'o(k) are equal to one for all k, except those that coincide with the vectors of reciprocal lattice 
^q. At these points So always drops to zero, while So drops to zero only for k = 0 and remains different 
from zero for finite k. 


The parity correlation function for a finite system is given by 

N 




( 164 ) 
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Figure 14: Critical temperature of the ideal Bose gas in a three-dimensional lattice. Solid line is an exact n umeri cal so lution 
of Eq. | |167[ l. Dashed lines are approximate analytical results for small and large fillings described by Eqs. ( |170[ l and jmj, 
respectively. 


provided that li ^ I 2 . If li = I 2 , the first term on the right-hand side of Eq. (1641 is replaced by one. 
Approaching the thermodynamic limit, (li,l 2 ) vanishes as 


A(-i)- (li ^ h) « -^(ni)exp(-4(hi)) 


(165) 


Eqs. (159), (162), (165) show that the leading finite-size corrections scale as ^ (hi)/L“. 

6.1.3. Critical temperature for the condensation 

For the ideal Bose gas, the total number of particles in the grand-canonical ensemble is given by 

1 


E 


kelBZ 


exp [(ek - m) /knT] - 1 ' 


(166) 


The critical temperature Tc is defined by the condition p = eo |171L 1188] . In the limit of infinite lattice, the 
summation in Eq. (166) can be replaced by an integral which leads to the following equation for the critical 
temperature: 


... N ^ 


t=i 


exp 


2J 

feeTc 


/ 2J . 


1 d 


where 


1 r 

Io{x) = — / exp {x cos ip) d(p 
27r 


(167) 


(168) 


is the modified Bessel function of the first kind. At large values of the argument, the asymptotics of Io{x) 
is given by 

Io{x) « exp{x)/V2Trx , (169) 

which implies that the series in Eq. (167) converges only for d > 3, provided that Tc is finite. Therefore, in an 
infinite homogeneous lattice the finite-temperature Bose-Einstein condensate exists only in three dimensions 
similarly to the case without any external potentials. 

Numerical solution of Eq. (167) presented in Fig. 14 show s that Tc grows with the filling factor (fii) = 
N/L^. For (hi) <C 1, Tc is small and Iq{ 2Jj/k-QTf) in Eq. ( |167[ ) can be well approximated by the asymptotic 
expression (169). This gives the result 

fcsTc 


J 


= 47r 


(^ 1 ) 

LC(3/2)J 


2/3 


C(3/2) = 2.612... 


(170) 


which can be also obtained from the well-known expression for the ideal Bose gas in a homogeneous contin¬ 
uous space |188| replacing the mass M by the effective mass (93) which follows from the dispersion relation 
of non-interacting particles in a lattice (140). 
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Figure 15: (color online) Energy eigenvalues (a) and the level spacing (b) for a single particle in a lat tice o f L = 200 sites in the 
presence of harmonic potential with j J = 0.005. Black dotted lines show the results given by Eq. ( |175[ | and red dashed lines 
are determined by Ei = which corresponds to the local energies of the harmonic confinement. (For interpretation of 

the references to colour in this figure legend, the reader is referred to the web version of this article.) 


In the opposite limit, (ni) ^ 1, is large and the sum in Eq. (1671 can be replaced by an integral. In 
this manner, we obtain 


k^Tc _ 2(ni) + 1 

~ ~ 73 


h 



[exp{-x)Io{x)f 


dx . 


(171) 


The integral I 3 was calculated in fact by G. N. Watson in 1939 nog and the result can be expressed in 
terms of the complete elliptic integral of the first kind 


as UnilM] 


K{x) 



(1 — x^ sin^ t) dt 


h 


f 18 + I2C2- I0C3 
\ 



(172) 


(173) 


which results in the numerical value I 3 = 0.505.... In the case of unit filling, (hi) 
k^Tc/J = 5.591 |204| . and for (hi) = 2 we get k^Td J = 9.69. 


1, Eq. (1671 gives 


6.1.4-- Harmonic trap 

In the presence of confining potential the translational invariance is broken. Since in the case of harmonic 
trap the problem is separable, we restrict ourselves to one dimension. In addition, we shall consider open 
boundary conditions. 

In the special case Vt = 0, the solution of the eigenvalue problem is given by 

i = (174) 

In the case of nonvanishing Vt the single-particle eigenstates become qualitatively different |205H208] . They 
can be separated into two parts and some transient regime between those. The lowest part of the energy 


Ei = —2 J cos 
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Figure 16: Eigenmodes of a single particle in a lattice of L = 200 sites in the presence of harmonic potential with VpIJ = 0.005. 


spectrum is well described by the textbook result for the harmonic oscillator in continuum (see Fig. 151. In 
terms of the parameters of the Bose-Hubbard model it reads 


E, = -2J + 2^/J¥T[i-]^ 


z = l,2,..., 


(175) 


and the agreement becomes better for smaller values of Vt / J■ The corresponding wavefunctions are localized 


in the middle of the trap, see Fig. 16 'a,b,c). 


If the energies grow, the level spacing is not constant anymore. For Ei > 2 J, the states become doubly 
degenerate which results in the rapid oscillations of the level spacing |206| . Fig. ^h). This part of the 
spectrum is determined mainly by the local energies of the trapping potential because the contribution of 
the hopping term is suppressed. The corresponding wavefunctions vanish at the trap center and oscillate 
in the regions between the classical turning points ic = io ^ \l{Ei + 2 J)/Vt and the turning points is = 
£o ± \/{Ei — 2J)/Vt associated with Bragg reflection [20511206] . With the increase of energy these regions 
move towards the trap edges and shrink which can be seen in Fig. 16 d,e,f). The Bragg reflection modifies 
significantly the density of states of a single atom even in the limit Vt —>■ 0 compared to the homogeneous 
case Vt = 0 with periodic boundary conditions: the square-root singularity in one dimension is replaced by 
a logarithmic one, and the logarithmic van Hove singularity in two dimensions disappears altogether 


6.2. The limit of vanishing tunneling 
6.2.1. Eigenstates 

In the absence of tunneling (J = 0), the eigenstates of the Hamiltonian are given by Eq. (641 and the 
corresponding eigenvalues 

^ {nn - 1) (176) 

1 

do not depend on K. The properties of the eigenstates appear to be different for commensurate and 
incommensurate fillings. 

We consider first the case of the commensurate Ailing, N = nL‘^, where n is a positive integer. The 
ground state has equal occupation numbers at each lattice site, that is nri = n and, therefore, the particle- 
number distribution is a Kronecker delta. Since any translation maps the state into itself, we have izr = 1 


in Eq. (641 which means that such a state exists only a,t K — 0. It has the energy 


pOO 
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This is a perfect insulator = 0) and all spatial correlation functions vanish for this state. 

The excited states are degenerate and form flat energy bands. The lowest band contains, at each value 
of K, — 1 degenerate eigenstates with the energies EU , T = 1,..., — 1. These states 

correspond to bosonic configurations with the same occupation numbers n at any site except two, one of 
which contains n — 1 bosons (hole excitation) and the other one n + 1 (particle excitation). They have 
vy = L for each spatial direction and, therefore, exist at any value of K. The highest band contains L'^ 
degenerate states (one state for each value of K) with all atoms sitting at the same lattice site. These states 
have the energy = UN{N — l)/2. As we will see later, finite hopping rate J lifts the degeneracy, the 
bands acquire finite widths and can even overlap if the tunneling parameter is large enough. If the Ailing is 
incommensurate, not only excited states but also the ground state are degenerate. 

In the case of large interaction but incommensurate Ailing N = nL‘^ + IV', N' < L‘^, the ground state is 
degenerate. These are states which are obtained from the state |'0nL'i) creating one additional 

particle on N' sites. Their energy is 


En = {L^ - (n - 1) + N'^ {n + l)n. (178) 

The degeneracy is lifted at least partially, if we switch on infinitesimally small hopping J. The energy of 
the first excited state will be infinitesimally close to that of the ground state, i.e., there will be no energy 
gap in the excitation spectrum. 


6.2.2. Finite temperature 

The grand-canonical partition function in the limit J = 0 has the form Z (/i) = (fj.) |209| , where 




n—0 


En - nn 
k^T 


U , 

En = - 1 ) 


(179) 


is the partition function of a single lattice site. From this we can deduce the on-site particle-number 
distribution 


p{ni = n) 


-J , X exp 


En - p.n \ 
k^T ) 


(180) 


At finite temperature, there is always a one-to-one correspondence between the chemical potential /i and 
the Ailing (hi) determined by the equation 


(«i) = X! ■ 

n—0 


(181) 


In general, this equation has to be solved numerically but in some special cases it allows also analytical 
solutions. 

If the temperature is low, k^T U, the particle-number distribution is strongly peaked near n = no, 
where no is an integer which is close to (hi), and the probabilities of the occupation numbers diflerent from 
no, no ± 1 are negligible. In this case, we obtain the following expression for the chemical potential 


exp 


(— 

UbT 


2 (1 -I- no - (hi)) 


(hi) - no-I-W ((hi) - no) +4 1 - ((hi) - no) 


(182) 


which generalizes the solution obtained in Ref. |210| for flllings (hi) close to one. The nonvanishing proba- 
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Figure 17: (color online) Probabilities p{ni = n) of the occupation numbers n = 0 (red), 1 (green), 2 (blue), 3 (magenta), 
4 (brown) for J = 0 and filling (ni) = 1 (a), (fii) = 0.5 (b) obtained according to Eq. (|180[ l with the aid of numerical solution 
of Eq. (|181|l for the chemical potential p,. Dashed lines are the low-temperature limit (|183|l with no = 1. 


bilities are given by 


p{n\ = no) = 

p{n\ = no - 1) = 


(1 +no - (ni))exp (^) 


(fc^) 
(1 + no - 

+ 2 exp 

(hi)) exp 

'U{no-l)' 

kuT 

'U{no-l)' 

k^T 

exp 

{ksT J 

+ 2 exp 

'U{no-l)' 

ksT 


(183) 


p(ni = no + 1) = 1-p(ni = no)-p(ni = no - 1) . 

If the filling is integer, (ni) = no, p, = U{nQ — 1/2) and these expressions simplify as 

1 , , exp(-^) 


p(ni = no) = 


1 + 2 exp (— 


( _’ 

V 2feBT J 


p{ni = no ± 1) = 


1 + 2 exp (— 


( _ IL_] 

V 2 ksT J 


(184) 


In the high-temperature limit, k^T ^ U, and in the case of arbitrary filling, the chemical potential is 
negative and has the form |209| 

p = -k^Tln ^ 

In this regime, the particle-number statistics is described by the geometric distribution 

(ni) 1 


p(ni = n) = 


1 + {hi) J 1 + (ni) 


n = 0,1 ,... 


(186) 


The complete temperature dependence of p{n\ = n) in the case of unit and half filling is shown in Fig. 17 
see also Ref. 


6.3. Scattering and bound states of two interacting atoms 

We consider two atoms in a one-dimensional lattice under periodic boundary conditions. We assume 


that L is odd and look for the eigenstates of the Hamiltonian (541 in the form of the superposition 


(L-l)/2 

\Kn) = Ciforln/fr) > 

r=o 


(187) 
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where the basis states (641 are explicitly given by 




VI 

1 


j=0 

L-1 


tk 


(188) 


L-l 


^ TTr'^i^j ^ ■ ■ ■ > 2 

i=o \ / r-i i-r-i 


where the index T in this particular system has a meaning of the interatomic distanc e an d we have introduced 
the notation \ni.. .ni) = ^£^1 \ni). The eigenvalue problem for the Hamiltonian (541 can be written down 
in the form 

(L-l)/2 

CKiir' = E 2 ^ckq.t ■ (189) 

r'=o 

The nonvanishing entries of the tridiagonal (L + l)/2 x (T + l)/2 matrix are given by |211| 


H^k^ = U, 

j^{L—l)/2,{L—l)/2 _ _ j 

Hr — O 

K" = = -JV2{1 + tk) , 


(190) 


(L+l)/2 

'K 


JL-l)/2 

'K 


H 


r,r+i 

K 


= -j{\ + tk) , r = i,... 


L- 3 
2 


The eigenvectors in Eq. (1891 satisfy the normalization condition 

= 1 . 


(L-l)/2 


E 

r=o 


Cicnr 


(191) 


In the case of even L, the summation in Eq. (187 \ is over T = 0,..., L/2. The basis state with T = L/2 can 
be rewritten in the form 


|nic,r=L/2) — E |10..^1^.^) 




Tk 


(192) 


which s hows that it exists only for even values of K. The form of matrix Hk remains almost the same as 


. . r r 

in Eq. (|190|). The only difference is that now all diagonal elements Hj^ with T > 0 vanish and the last 

00 


off-diagonal matrix element 


H 


L/2-l,L/2 

K 


= H 


'Jq,2m 


(193) 


m— — oo 


for even q is the same as 


The eigenvalue problem (1891 was solved in Refs. |211H213] analytically and numerically for negative U 
but the result can be easily generalized to arbitrary U. Later it was discussed in the context of ultracold 
atoms The energy spectrum consists of the scattering states of a pair of asymptotically free 

particles and the bound state. 

The energies of the scattering states are given by Eq. (1531. In the limit L —>■ 00 , they form a con¬ 
tinuous band shown in Fig. 18 with the boundaries Ek,o = —Qk and Ex^-K/a = <1k- The corresponding 
eigenstates (1871 are describeddDy the coefficients 


Cif.fc.r = CK,k,o 


72 


cos(fcoT) 


U sin (fcoT) 


qk sin (fca) 


Ka^ 

exp 1 


(194) 
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Figure 18: Energy spectrum of two interacting atoms. Shaded region is the scattering continuum described by Eq. | |153[ l. 
Solid line below the scattering continuum and the dashed line above it are the energies of the bound state ( |196| l with U/ J = —5 
and V! J = +5, respectively. 



Kajir 

Figure 19: Mean distance between two atoms in the bound state ||197[l with U/J = ±5. 


and the normalization has been discussed in Ref. |220| . In the limit L —>■ oo, the bound state is given by 


CKO = 


CRT = 


'1^ 

l + bl^ 


bK = 


U-S 


K 


QK 


(195) 


V^CRob^ exp 


Ka^ 


r = 1,2,..., c» 


and has the energy 


E. 


KQ, 


= £r = U\ 1 


(f). 


(196) 


which is also shown in Fig. [T^ 

In the absence of interactions {U = 0), bR = ±1 and the normalization condition (1911 cannot be fulfilled 
which means that the bound state does not exist in this case. As soon as the interactions are present {U ^ 0), 

I < 1 which guarantees the normalization. 

In the case of attractive interactions {U < 0), bR > 0 and all the coefficients crt are positive. The 
wavefunction does not have nodes since we are dealing with the ground state of the system. In the case of 
repulsive interactions {U > 0), bR < 0 and the sign of the coefficients Cry alternates. The wavefunction has 
infinitely many nodes reflecting the fact that this is a highly excited state. 

The distance between the atoms u; is a random variable which take s the values ic = T, with the proba¬ 


bilities |ciyr| ■ The mean interatomic distance in the bound state (1951 is given by 


2/)2 

(wr) = ^ 


l-b% 


(197) 


which is shown in Fig. 19 (wr) takes its maximal value at AT = 0 but vanishes at Ka = tt. 

Momentum distribution appears to be drastically different for different types of the eigenstates. For the 
scattering states, it has two sharp peaks corresponding to the momenta of the atoms. The bound states 
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Figure 20: Quasi-momentum distribution of two interacting atoms in the bound state with K = 0 and VjJ = —5 (solid), 
+5 (dashed). 


are characterized by broad momentum distributions shown in Fig. |20[ In the case of attractive interactions, 
the momentum distribution takes its maximal value at A: = 0, while for repulsive interactions, the maxima 
appear k = -in:ja. With the increase of the interaction strength |(7|, the momentum distribution of the 
bound states becomes broader |221| . 

The characteristic feature of the bound state is that \ckq\ > |cifr| > T = 1,...,(L — l)/2, i.e., the 
probability of finding two atoms on the same lattice site is higher than all the other ones. In the case of 
attractive interactions {U < 0), this sort of localization corresponds to the soliton solution of the discrete 
nonlinear Schrodinger equation and, therefore, the discrete level can be called a “soliton band" mu. 

The existence of the bound state in the case of repulsive interaction {U > 0) is quite unusual from 
the point of view of classical physics. This is a purely quantum phenomenon which occurs in a structured 
environment produced by the periodic potential in the absence of dissipation. Repulsively bound pairs of 
®^Rb atoms were created in the experiment m using a magnetic-field sweep across a Feshbach resonance 
and demonstrated by measurements of the momentum distribution and the binding energy. 

It is interesting to compare the solutions obtained for identical bosons with those for distinguishable 
particles. In this case the energy spectrum remains the same, although the states are formally different. 
They are linear combinations 


L-l 


\KQ) = ^ c/cnrlnicr) 


(198) 


of the basis states 


Injco) = 


In^r/ = 


r=o 






L-l 

r -1 L-r-i 


r = i,...,L-i 


where the indices 1 and 2 label the two atoms. The bound state in the limit of large L is described by the 
coefficients 


CKF = CKob^K exp ( i^F 


Ka, 

Y 


ck,l-t = c*j^T > F = 0 , 1 ,..., Ll2 , 


(199) 


where ckg is given by Eq. (1951. In spite of this difference, the distribution of the interparticle distances, 
which is given now by |cico| for F = 0 and |cicr|^ + \ck,l-t'^ for F = 1 ,..., T/2 , remains the same as in 
the case of the indistinguishable bosons. 


6-4-. Hard-core bosons 

In the limit of infinite repulsion {U = oo), the occupation numbers of the individual lattice sites are 
restricted by 0 and 1. Formally, this leads to the constrains for the bosonic operators = 0 and 
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|aj, dj| = 1 and also implies that the total number of atoms N cannot be larger than the total number of 

the lattice sites L‘^. The state with N = L‘^ is trivially an insulator as no hopping on individual sites can 
take place. Non-trivial physics is possible only for N < L‘^. 

In this limit, bosonic operators can be mapped into spin-1/2 operators by means of the Holstein-Primakoff 
transformation mg 


at = al 


1 — a 


i“i > 


= \ 1 — d, d 


i“i “1 


0-1 


= 2d|d[ — 1 


( 200 ) 


where a^ = (Tj“ ± ia\ are raising and lowering operators and d/, a = x,y,z, are Pauli matrices. In what 
follows we assume that the products of the bosonic operators dj, dj are arranged in the normal order, 
that is, c reation operators are placed to the left of the annihilation operators. Then the square roots in 
Eq. (200) become identities and the bosonic operators can be directly replaced by |1851 1186] . The 
transformation (200) maps the Bose-Hubbard Hamiltonian into the XY spin-1/2 Hamiltonian: 


Hxy = — 2 J 


d 

EE 

1 i/=i 


cr, cr; 


1 ‘^ 1 +e. 


+ 


1 "n-e, 


J + 


E.A 

1 


1 


( 201 ) 


6.4-.1. Bose-Fermi mapping in one-dimensional lattices 

The limit of infinite repulsion of bosons in a one-dimensional lattice is called Tonks-Girardeau regime. 
It is exactly solvable via the Jordan-Wigner transformation 

di = exp ( i-K ^ c\,c^, I (l - 2 cj,c^,^ ci , (202) 

V e'<e J e'<e ^ 


where q and c\ are the fermionic annihilation and creation operators. Under this transformation, the 
hopping term of the Bose-Hubbard Hamiltonian as well as the local particle-number operators djd^ remain 
invariant and the Hamiltonian (54) in the presence of an arbitrary external potential ei is mapped onto the 
one of noninteracting fermions: 


Hbh — Hf — J ^ ^ 


eec\ci . 


(203) 


i=i 


Periodic boundary conditions for bosons are equivalent to the requirement 

CL-i-i = exp Cl , 


(204) 


which implies periodic boundary conditions for fermions if the number of particles N is odd, otherwise one 


should use antiperiodic boundary conditions in the Hamiltonian (2031. This feature requires some care in the 


studies of the hard-core bosons with periodic boundary conditions in the grand-canonical ensemble because 
the states with even and odd particle numbers should be treated in general separately |225[ 1226] . In the 
case of open boundary conditions for bosons, the boundary conditions for the equivalent fermionic system 
remain the same. Comparison with exact numerical solutions for the one-dimensional Bose-Hubbard model 
for arbitrary interaction U shows that the fermionization approach gives quantitatively correct results for 
U/J> 200 [2271. 


The Y-particle eigenstates \'ijjp) of the Hamiltonian (203) can be constructed as products of the single¬ 
particle eigenstates 


«) = 4 | 0 ) = 'yipaec] | 0 ) 


(205) 
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with the eigenenergies Sa, a = 1,... ,L. If the energies are labeled in the ascending order, i.e., ei < £2 < 
• • • < El, the ground state is given by = na=i lo^)- 

The superfluid stiffness is deflned by Eq. ( [M| ) . In the limit 0 —>■ 0 and at T = 0 it is given by Eq. (96 1 . 
For hard-core bosons in ID, the latter takes the following form m 


N 




i: E 


N 




1 


i—1a—1 


N 


a=N+l P=1 


,i+i ~ ‘Pa,e+i‘f/3e) 


(206) 


where ipa,L+i = (- 1 )^+Vai. 


From Eq. (202), it follows that the bosonic Lx L one-body density matrix with the entries {a\af^,) cannot 


be simply identified with the fermionic one (cjc^,). Although for i' = {of^a^,) coincide with {c\c^i), in 

general they are not the same. However, they are related to each other and for I' > I the quantities (ajo,^,) 
can be worked out as determinants of the {£! — €) x {£' — £) matrices ^ [H 12241122511228L1229] : 


{a\ag,) = 2^ ^ ^detG^^’^^ 


(207) 


where the entries of ^ are given by 


Gff ^ , f, j = 1, ...,(£'- ^) . 


(208) 


At zero temperature, fermionic one-body density matrix (cjc^,) can be calculated using the solution of the 
single-particle eigenvalue problem as 


N 

a=l 


If the temperature is finite and for open boundary conditions, Eq. (2091 is generalized by 

L 


( 4 g ') = E siev’d'fFD{ea) 


(209) 


( 210 ) 


where /fd is the Fermi-Dirac distribution function: 

/fd(£) = 


exp [(£ - fj.) //cbT] -I- 1 


( 211 ) 


The chemical potential /r is fixed by the requirement to have the desired number of particles N = 

Alternatively, the bosonic one-body density matrix at T = 0 can be represented in the form [28011281| 


(aja,,) =det [P^(f)P(£)] 


( 212 ) 


where the L x {N + 1) matrix P{£) is determined as 


1 

f 

for 

i = 1,.. 

.,£-1, 

a = 1,...,A^; 




for 

i = £, . . 

.,L, 

a = 1,...,A^; 

(213) 

1 

[ he. 

for 

i = 1,.. 

.,L, 

a = N + 1. 



Eq. (2121 can be also rewritten in the form |282| 


N 


{a\a^,) =dei[A{£J')] y] 

Q ,^—1 


(214) 
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where the entries of the N x N matrix A{£,£') are given by 


e'-i 


= l<£! . 


( 215 ) 




The advantage of Eqs. (2141, (2151 is that the numerical calculation of (aja^,) can be done efficiently thanks 
to the recurrence relation 

A^p[£,£' + 1) = A^0{£,e) - . (216) 

It is interesting to note that Eqs. (214), ( 215[ ) are the discrete version of the ones for the Tonks-Girardeau 
gas in the continuum [233]. 

For finite temperature and in the grand-canonical ensemble with open boundary conditions, the bosonic 
one-body density matrix for £' can be also obtained as [226] 


(aja^,) = - {det [/ -f (/ -k A)0{(!)UR U^O{£)] - det [/ -k 0{£')UR t/1'C>(£)]} , 


(217) 


where 


L r 


Z = det(/-ki?) = 


a=l 


1 + exp - 


£a- fJ- 

k^T 


(218) 


is the partition function, I is the identity matrix, U is the unitary L x L matrix containing the eigenvectors 
of H-p for = 1 in its columns {Uia = 0 (£) is a diagonal matrix with the first {I — 1) elements equal 

to —1 and the others equal to -kl. In addition, we defined R = exp [—(e — iJ,I)/k-QT\ with e being a diagonal 
matrix of the eigenenergies £„. The diagonal elements of the one-body density matrix are the same as for 
noninteracting fermions and can be easily calculated according to Eq. ( |210| ) with £! = £, which can be also 
rewritten as [333] 


(aja,) = (cjc,)=[[/(/ + 

After simple transformations, Eq. ( plTt can be rewritten in the form analogous to (|214|) [232] 




(219) 


L 

CK,^—1 


( 220 ) 


where the Lx L matrix B{£,£') is defined by 




( 221 ) 


and A is the same as in Eq. (2211 but extended to a, /3 = 1,..., L. 

The particle-number correlation function of hard-core bosons is formally the same as for non-interacting 
fermions: 

Fni^A') = {nf)5u'- (c|c^,) . 


( 222 ) 


At finite T and under open boundary conditions, it has the form 

L L 

FniiJ') = Y ‘PatfFDiSa) Y “ fpoiefi)] 


(223) 


9=1 


which allows to express the dynamic structure factor in terms of the single-particle eigenmodes as [333] 


S{k,uj) = Y 

a ,0 


E 


‘Pae^pee 


Akai 


/FD(ea) [1 — fFDisp)] <5 ( W — 


£ (5 £c 


(224) 
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The parity operator (1461 has a form = 1 — 2cjc^ and, therefore, the parity correlation is (^i,^ 2 ) = 

4F„(£i,£2)- 

Calculation of the noise correlations reduces to the calculation of the one-body density matrix and of the 
four-point correlation function The former has been already discussed above and the latter 

at r = 0 can be computed as |186| 


= det [P'''(4,^ 2 )] , 
where the L x {N + 2) matrix P{£i,£ 2 ) is given by 


1 

r -pui2) 

for 

i = l,.. 


cx — 1 ,..., + 1; 

PU£i.^2 ) = 

P^a{h) 

for 

i = £i,. 

..,L, 

a = 1,..., 

1 


for 

1 = 1,., 

,.,L, 

a = A^-h2; 


(225) 


(226) 


and Pia{£) are determined by Eq. (2131. A different approach to the computation of noise correlations based 
on Wick theorem was developed in Ref. |185| . However, it appears to be less efficient. 

Finally, we would like to note that the time evolution of an arbitrary initial state 


N L 


iv’f(o)) = n XI 

t=l 


(227) 


can be easily calculated P55] 





N L 

=n E 4 |0) , 

OL—l E—1 


(228) 


which remains to be a product of single-particle states similarly to the ground state IV'f)- difference 

is that the coefficients ipai become time-dependent: 


L L 

li = l p=l 



(229) 


Therefore, the time evolution of the observables is described by the same equations as for the expectation 
values in the ground state (at T = 0) but with ipat replaced by ipatit)- 


6.4-2. Homogeneous lattiee 

Now we apply the general formalism described in the previous section to the homogeneous lattices (e^ = 


0). In this case, the solution of the single-particle problem has the form (1491, (1501 with kaO = 2qa'!T/L, 


go, = 0 ,..., L — 1, for periodic boundary conditions, and = (2qa + 1 )ttIL for antiperiodic. 
The ground-state energy has the same form for even and odd N: 


A™ = - 2 J 


sin(7r(h£)) N 

■ , iT\ ’ = ~T ■ 

sm(7r/L) L 


(230) 


In the thermodynamic limit, this leads to the following expression for the chemical potential: 


p =—2Jcos {tt{ n^)) (231) 

and, therefore, the compressibility is given by 

K = [27r Jsin (7r(h^))] ^ . (232) 

Note that fi decreases with J for 0 < (h^) <1/2 but increases for 1/2 < (n^) < 1. 
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The energy spectrum of the iV-particle system in the thermodynamic limit is continuous and has no gaps 
with the minimal and maximal energies equal to ±2 JL sin {TT{ni}) /tt. The lowest excitation mode at small 
momentum has a linear dispersion relation characterized by the sound velocity |236| 


Cs = 2a— sin (7r(nf)) 


(233) 


which coincides with the Fermi velocity 


1 dek 


k—kp 


where the Fermi momentum fcp is defined as 


kp = k„^Ni2 = - {hi) . 
' a 


(234) 


(235) 


The superfluid stiffness is described by Eq. (2061. In the case of a homogeneous lattice, the second term 
vanishes and we obtain 


fN = 


sin {TTihi)) 


^ iVsin(7r/T) 

In the thermodynamic limit, this gives the well-known result (see, e.g.. Ref. |237p 


(236) 


f = 

J OO 


sin(7r(hf)) 

Tr{hi) 


(237) 


From Eqs. (2321, (2331 and (2371, one can easily see that the superfluid stiffness is related to the sound 

M*Csa 


velocity as |238] 


r = 

J OO 


Th{hi) 


(238) 


and the hydrodynamic relation (95) is also fulfilled. 


In the ground state, the onsite particle-number statistics is described by a binary distribution 


p{ni = n) = (1 - {hi)) 6n,o + {hi)Sn,i ■ 


(239) 


This leads to the result for the particle-number fluctuations 

<Jni = V {hi) (1 - {hi)) 


(240) 


which follows also from Eq. (222). 


The fermionic one-body density matrix defined by Eq. (209) is given by the explicit expression 


_ sin [K{hl){^ - ^’)\ 
~ Lsin[f(£-f)] 


(241) 


which appears to be the same for even and odd N. Therefore, matrices ^ needed for the calculation of 
the bosonic one-body density matrix in Eq. (2071 acquire a Toeplitz form 




^ 9i 

90 

9i 


9e'-e-2\ 

92 

9i 

90 



93 

92 



9i 




91 

9o 

\9£'-i 


93 

92 

91 / 


(242) 
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Figure 21: One-body density matrix of hard core bosons in a one-dimensional homogeneous lattice with periodic boundary 
conditions in the thermodynamic limit for N/L = 0.5, 0.3, 0.2, 0.1 (from top to the bottom). Red dashed line shows the 
asymptotics at large distances in the case of N/L = 0.5 according to Eq. ( |247| l. (For interpretation of the references to colour 
in this figure legend, the reader is referred to the web version of this article.) 


where 


sin(7r(n^)n) 1 
9n = ~ : 7^ ^ 


(243) 


In the case of half filling, (n^)=l/2, the coefficients Qn with even n vanish and the bosonic one-body 
density matrix is given by 


{a\a^,) = 


— 7?^ 

^ “ 2 “ 

2 R e'-e-1 R t'-t+i 


for even I' — 
for odd — t, 


(244) 


where 


Rq=\- 


q 9-1 

n 

fc=i 


sin^ {2k'K/L) 


sin [{2k -I- 1 ) tt/L] sin [{2k — 1 ) tt/L] 


q—k 


In the thermodynamic limit, R„ simplifies as 


f-Vii 

{2kf 


{2k + l){2 k-l) 


q—k 


(245) 


(246) 


At large distances \f — i\, the bosonic one-body density matrix has the following asymptotics |23911240] : 


{a\a^,) 


C 




1 - 


(i)K'-d 


C = 0.294176... 


(247) 


The dependence of (a)a^,) on the distance ]£■' — P\ in the case of half filling described by Eqs. (2441, (2461 
is shown in Fig. 21 Surprisingly, the asymptotic formula (2471 is in perfect agreement with the exact result 
already at small |£' — ^\. For other fillings, the correlation function {dla^,) is smaller and has a different 
behavior at small distances. However, at large distances it shows the same asymptotics {d\d^,) ~ \£' — 
which is a manifestation of quasi-long-range order. This implies that the quasi-momentum distribution has 
a singularity at fc —J. 0 |241fi243] . The largest eigenvalue of the one-body density matrix which 

corresponds to the zero quasi-momentum state scales as y/N |231| . 
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kaj'K 


Figure 22: (color online) Static structure factor of hard core bosons in a one-dimensional homogeneous lattice of L = 100 sites 
with periodic boundary conditions. N = 40, 50,60, 70, 80, 9 0.100 f rom top to the bottom. The data obtained at discrete p oints 
kq are connected by straight lines. So{kq) is given by Eqs. (248 i, ( 249 i. So{kq) is obtained from So{kq) according to Eq. (129 I 
using Goik) for Vo = 10Er. 


Any finite temperature has a strong infiuence on the properties of the one-body density matrix |226| . 
The long-distance asymptotics becomes exponential, ~ exp (—|^i — ^2|/0i the quasi-long-range 

order present in the ground state is destroyed. The correlation length ^ decreases with T and for small 
temperatures decays as ^ ~ 1/T. It has also strong dependence on the filling. 

The static structure factor 5'o(fc) is a periodic function of k with the p eriod e qual to the vector of the 
reciprocal lattice kq^L = It can be easily obtained at T = 0 from Eqs. (1301, (|209|), (2221. Within one 
period and for N = 1,..., L/2 it is given by 


Soikq) = 


kq 

2 /cp 

1 

kL—kq 

2 kp 


7 kq 0, . . . , 2/tp , 

, q = 2kp ,... - 2fcF , 

, kq — kij 2kp ,..., kij , 


(248) 


where kq = j^q and fcp is the Fermi momentum. In this case, > Akp. For N = ^ + 1,... ,L we have 


instead 


Soikq) = 


^ , kq = 0,... ,kL - 2kp , 

, kq — k]^ 2kp ^..., 2kp 


2 Alp 

_ ) kL- 2 kp 
2kp 

kL — kq 

2fcp 


(249) 


, kq — 2 kp ,..., kij , 


with 2 A:f < fci < 4 A:f- For N = 1,..., L — 1, Soikq) is a linear function of kq in finite intervals in the vicinity 
of fc = 0 and k = 2'n/a with the slope l/(2fcF) in complete agreement with Eqs. (145), (232), ( |233 1. In the 
remaining interval near fc = ir/a, Soikq) takes a constant value which is one for iV = 1,..., L/2 and reduces 
from iL — 2)/(L + 2) to zero for = ^ -|- 1,..., L. As it follows from Eq. (249), Soikq) vanishes ioi N = L 


and the s tatic structure factor in continuum takes the form S'o(fc) = 1 — (j'g(fc), where Goik) is defined by 


Eq. (132). The comparison of S'o(fc) and Soik) for different values of N is given in Fig. 


22 


Noise correlations of hard-core bosons in homogeneous lattices possess the following characteristic fea¬ 
tures |18511186| . (i) Very large peaks appear at fci = ^2 = 0 as well as at integer multiples of the reciprocal 
lattice vector as a result of the quasi-condensation and the underlying order induced by a periodic lattice. 
The heights of the peaks depend on the filling, (ii) A line of maxima exists for ki = ^2 due to the bunching 
typical for bosonic systems, (iii) There are dips immersed in a negative background along the lines fci = 0 
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Figure 23: (color online) Spatial distribution of the mean occupation numbers of hard-core bosons in a one-dimensional 

lattice and in the presence of harmonic confinement with the parameter J/Vr = 200 at T = 0. The total number of particles 
N = 1, 2, 3,4, 5,10, 20, 30, 40, 50, 60, 70, 80, 90,100 from bottom to the top. (a) Exact calculations in a lattice of L = 200 sites 
(the size of the lattice does not affect the results), (b) Local-density approximation ||251^. 


and ^2 = 0 which are related to quantum depletion, (iv) The correlation function for ki = k 2 = 0 scales 
linearly with the system size. 

6.4-3. Harmonic trap 

We consider the effects of the harmonic trapping potential described by the local terms 


Q = Vt (^ — ^oY 


r T = -- 


(250) 


in the Hamiltonian (2031, where £o denotes the center of the trap (which is not necessarily an integer). 
Spatial distributions of the mean occupation numbers (nY) for different total particle numbers N at T — 0 
determined by Eq. (2091 with £' = £ are shown in Fig.[^ With the increase of N, (h^)’s become larger and 
the size of the atomic sample grows. If the particle number exceeds the value N « 2.Qd>^JJ/V-£, a plateau 
with {ni) = 1 appears at the trap center |2061124411245] . In the example shown in Fig. 23 this happens for 
TV > 38. Since the local chemical potential p-i = p, — varies from site to site but the occupation numbers 


(■hi) within the plateau do not, the local compressibility (89) vanishes. 


The density profiles can be also obtained within the local density approximation jH 124611247] replacing 
the chemical potential p in Eq. (231) by = /r — e^. This leads to the following expression 


{nY) = \ iarccos(-^) 


; P^ 2T, 

> \k‘i\ ^ 

; /4£ ^ 2T, 


(251) 


where the global chemical potential p is determined by the total number of particles as ^■ ^ig- 23 

shows that the local-density approximation works very well for N > 10, although it does not describe fine 
details of the exact density profiles. 

One-body density matrix as well as natural orbitals in the presence of a trapping potential were studied 
in Refs. |2291I230] . It was shown that the power-law decay (aja^,) ~ \£ — £'\~^/'^ is preserved for intermediate 
distances but at larger distances the correlations drop much faster |230| which removes the singularity in the 
quasi-momentum distribution for /c —>■ 0 |229| . The harmonic trap sets a momentum scale px = '/MHtJY 
below which the momentum distribution is flattened due to suppression of the long-range correlations [Ml 
1227] . The presence of the plateau in the density profile is accompanied by the splitting of the natural orbital 
of the one-body density matrix with the largest eigenvalue into two parts localized in the side regions with 
nonvanishing compressibility piTj . 

In the experiment of Ref. |5] , a two-dimensional array of independent one-dimensional chains with the 
filling factor smaller than one has been created. Due to the harmonic confinement, the number of atoms in 
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a chain labeled by a pair of indices (i, j) is given by 


N,j = N> 


00 


1 5-/V ,, ON 


3/2 


where N is the total number of atoms and Nqq is the number of atoms in the central chain. The probability 
of having a chain with N' atoms reads 


P{N') 


2 


N' < Noo . 


The only parameter which determines this distribution was estimated to be Nqq = 18. The quasi-momentum 
distribution calculated for this setup using the Bose-Fermi mapping at finite temperature is in excellent 
agreement with the experimental data 1248] . 


6.4-4- Extended fermionization 

As it was mentioned above, standard Bose-Fermi mapping is not valid for N > L. Nevertheless, this 
formalism can be used within the framework of the extended fermionization model |1891I249] . The idea is to 
divide the whole system into two subsystems: (i) ub bosons siting at each lattice site, and (ii) N' = N — ubL 
excess bosons, where tt-b is such that N' < L. Then the excess particles can be treated as ordinary hard-core 
bosons with the effective tunneling parameter J’ = J{nB + !)■ This approach was employed to study the 
dynamic structure factor of hard-core bosons in the homogeneous one-dimensional lattices as well as in the 
presence of a harmonic trap |189| as well as damping of dipole oscillations |249| . In the homogeneous system, 
it leads to straightforward modifications of the results of section 103 For instance, the expression for the 
superfluid stiffness becomes 


riB + 1 sin {ttN'/L) 
N sin (n/L) 


(252) 


7. Perturbation theory in the limit of strong interaction 

If the tunneling parameter J is small compared to the interaction energy U, perturbative solution of the 
Bose-Hubbard model can be obtained in powers of J/U, which is called strong-coupling expansion |5I1 1^ . 
This method works very well for J/U below the critical point {J/U)c, which is less than one in all dimensions, 
where all physical quantities are analytical functions. It was used in the studies of spinless bosons with 
local IMS and nearest-neighbor interactions |26()l I2blj . two-species bosons |262] . and spin-1 

bosons |268L 1264] in different types of regular lattices like hypercubic isotropic |51L [52l [Ml 125511260H26311265] 
and anisotropic |253| . superlattices |25411255] . two-dimensional triangular and kagome lattices |256H258] . 
and in the presence of artificial gauge fields |25911265] . as well as in disordered lattices |521 1571 1266H268| . 
The strong-coupling expansion allows to avoid finite-size effects and shows excellent agreement with exact 
numerical data. 

However, this method can be applied only for commensurate fillings N = nL‘^, where n is an integer, 
or close to it for N = nL‘^ ± 1 due to the degeneracies of the zeroth-order eigenstates which may become 
intractable in the case of arbitrary filling. Analytical calculations for arbitrary d and n can be done in 
practice only for few lowest orders. Higher-order studies require symbolic calculations on a computer. This 
has been done for d = 1 and n = I up to the I4th order in Ref. m for the ground-state energy, variance of 
the occupation numbers and two-point correlation functions. Recently, symbolic calculations were done for 
d = 1 and n = 1, 2,3 up to the 16th order and the results were reported for the energy, two-point correlation 
functions and the first five moments of the on-site particle-number distribution [269] . Numerical calculations 
along these lines were also performed for d = 1, n = 1,2, and for d = 2, n = 1 up to the 13th order |2511I252] . 
In Refs. |791 I501[270| . the method was developed such that it became possible to do calculations in principle 
for arbitrary d and n. With the aid of the scaling theory it is possible to extrapolate to the infinite order of 
J/U |52l l54l r251| . Here we present analytical results for one-dimensional and isotropic hypercubic lattices 
of arbitrary dimensions. 
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7.1. Ground state in the case of commensurate filling 

The ground state in the case of commensurate filling {N = nL'^) is denoted as As it was 

discussed in Sec. |6.2[ in the zeroth-order of the perturbation theory, this state is not degenerate and has 


the energy (1771. Therefore, in order to study the ground-state properties of the Mott-insulator, one has to 
employ non-degenerate perturbation theory. For the ground-state energy per lattice site at arbitrary d and 
n up to the 4th order in J/U, this gives |271| 


_n{n-l) ( J 


ULd 


- ( A j Zn{n + 1) - ( ^ 


U 


U 


^n{n_p) _ 34 + ( 75 ^ _ 157)n(n -f 1 )] , (253) 


where Z = 2d is the coordination number. Symbolic perturbative expansion on a computer for d = 1, n = 1 
allowed to obtain the result up to the 14th order m- 


eT 

AUL 


J 
U 

4902596 /J 


6561 


jy 68 fJ . 

u) ’ 

12 


U 


1267 / J 
U 

14 


81 


44171 / J 


1458 \U 


10 


8020902135607 / J 
2645395200 [ U 


For d = 1, n = 1, Eq. (2531 reproduces the first two terms in Eq. (254). 


The probabilities to have ni particles on a lattice site are given by |271| 


(254) 


p{ni = n - 1) = ( ^ ) n{n + 1)Z 


J \ n{n + 1)Z 


U 


18 


p{n\ = n -I- 1 ) = 


J ^ n{n + 1)Z 
U 


18 


p(ni = n - 2) = ( — 


[8AZ - 156 -h n(334Z - 703) + ¥(326Z - 695)] 


n(n + 1)Z 


[76Z - 148 -h n{318Z - 687) + n^(326Z - 695)] 


i(n^ — 1)Z 




p{ni = n + 2) = 


n{n + l)(n -I- 2)Z 




and the probability p{n\ = n) can be obtained from the normalization condition 

n+2 

^ p{ni) = 1 . 

n\—n—2 


(255) 

(256) 

(257) 

(258) 

(259) 


Probabilities of other occupation numbers vanish in this order of the perturbation theory. The probabilities 
p{ni) given by Eqs. (255 l-(|2^l satisfy the relation 


p{n\ = n -I- 1 ) — p(n\ = n — 1) + 2 [p{n\ = n + 2) — p{n\ = n — 2)] =0 


(260) 


that follows from the obvious condition (rii) = n. Second-order terms in Eqs. (253), (255), (256) were 
obtained in Refs. |51[I5^[272] . 

Elements of the one-body density matrix ^ 0 ( 11 , 12 ) = (aqOb) third order in J/U have the 
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form [51] 


Fais = 1) = ^2n{n + 1) + [16Z - 34 + (76Z - 157)n(n + 1)] , (261) 

F,(s = 2) = 3n(n + l)(2n + l) , (262) 

Fjs = V2) = 2Fais = 2) , 

Fa(s = 3)=l^0 4n(n + l)(5n2 + 5n + l) , (263) 

Fa(s = V5)= 3Fa(s = 3) , Fa(s = V3)= 6Fa(s = 3) . 


Note that integer distances, s = 1, 2,..., are possible in lattices of any dimensio n d, w hile i rrational distances 
s = v^, and s = \/3 exist only for d > 2 and d > 3, respectively. Eqs. (253)-(261l can be tested for 
self-consistency using the identities 


ULd 

dU 


-:Lzf,{s = i) + 

(h 2 ) - (hi) 


{n'f) - (hi) 


(264) 


that follow from the Hamiltonian. 


With the aid of Eqs. (2611-(2631, one can deduce perturbative results for the quasi-momentum distribu¬ 
tion in the ground state: 


P(k) 


1 

N 


d 

n 4 - ^ 6 s n COS{S!yk,ya) , 

s ly—l 


(265) 


where bg = bsj^...s^ are invariant under any permutation of indices. In one dimension, the nonvanishing 


coefficients in the third order of the 

perturbation theory are given by 



hs = 

:2P,(s), s = l,2,3. 

(266) 

In two dimensions, we have 

bso = 

2Fa{s) , s = l,2,3, 

(267) 


bn = 

8Fa{2) , 6 i 2 = 12F,(3) , 


and in three dimensions: 




bsoo = 

2 Pa(s) , 

s = l,2,3, 

(268) 

&110 = 

8 Pa(2) , 

6120 = 12Fa(3) , 6111 =48Pa(3) , 



with all possible permutations of the indices. 


For the particle-number correlation function (1311 one finds [271] 


F„(l) = 2nin + l)(^^^ [64Z - 190 + n(n + 1)(448Z - 1069)] , 

/ 7 \ ^ 2n 

Fn{‘2) = ~ ( (11 + 43n-|-64n^-f 32n^) , 

F„(V2) = y(20-f 79n-4ll8n^-f 59n3). 


(269) 

(270) 

(271) 
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The consistency of Eqs. (269l-(271) with ( 253| ), (2611 can be also tested using the fact that the fluctuations 
of the total number of particles vanish: {N'^) — N'^ = 0. In a translationally invariant system this leads to 
the condition |269| 

L-l L-l 

^ ^ ((hohi) - (no)(hi)) = 0 . (272) 

^ 1=0 ld—0 


Considering only the terms up to the 4th order of the perturbation theory, we get 


n{n-l) J Z 


F„(l)+F„(2) + ^^F„(v^) 


= 0 , 


(273) 


where we have taken into account Eq. (2641. Eqs. (2531, (261), (269)-(271l satisfy indeed Eq. (273). 


The parity correlation is given by Em 

E(_i)r.(l) = ^ —^ 8n(n + l)+^—^ 


4n(n + 1) 


[16Z-34 + n(n +1)(40Z-157)] , (274) 

E(_i)„(2) = ^(7 + 29n + 44n2 + 22n3) , (275) 


F(_i).(y2) = ^(16 + 83n + 134n2+67n3) . 


(276) 


The results for the parity correlation function up to the second order in J/U at finite temperature can be 
found in Ref. UMl. A multi-point generalization of the two-point parity correlation function E’(_]^)n(R, I 2 ), a 
so-called string correlator defined by Eq. (2941, was calculated at zero temperature in the second order of the 
strong-coupling expansion for lattices of arbitrary dimensions and in the fourth order in one dimension P75] . 


1.2. Lowest excited states 

At each value of the discrete total momentum fiK, the first excitation band consists of — 1 states. 
In the limit J = 0, the states are degenerate and correspond to bosonic configurations with the same 
occupation numbers n at any site except two, one of which contains n — 1 bosons and the other one n -|- 1. 
In one dimension, they are explicitly given by 

= |nicn) = ^ |n-H, rt, .^. ,n ,n- l, n, .^. ,ri .) , (277) 


where LI = 1,, L— 1. In the first order of the perturbation theory, the degeneracy of these states at a given 
value of K is completely lifted. The energies of the excited states for one-dimensional lattices were calculated 
in the first order of JjU for unit filling n = 1 |25i)| and arbitrary filling |274| . By doing calculations for 
arbitrary n up to the second order in J/U, we obtain 


t^KQ. 

^nL 

U 


J 




nL 


J 


u 


u 


n 


-I- 1 - 2— cos ( TT— ] \lf + 4n(n -I- 1) cos^ ( — 


Ka 


(278) 


-!-( — ) 5n‘‘ -I- 6n -I- 2 — 4n(n -|- 1) cos ( 27r — 


n \ 2n{n -I- 1) -f (2n^ -|- 2n -|- l) cos Ka 


1 -I- 4n(n -I- 1) cos^{Ka/2) 


cos Ka 


where Lt = 1,..., L — 1. Eq. (278) has been derived for sufficiently large one-dimensional lattices {Z = 2) 
and can be easily generalized to arbitrary dimensions. The first-order term in Eq. (278) predicts symmetric 
form of the excitation band with respect to E = E^^ -\- U. The second-order term describes the asymmetry 
of the band. 
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The lowest excited state is labeled by itT = 0, = 1. From Eq. (2781 it follows that in the thermodynamic 


limit the gap separating it from the ground state (neutral gap) is given by 


U 


K\ - Kl 

u 


= l-2^(2n + l)+ (n" + 2n + 2) . 


(279) 


7.3. Particle-hole exeitations 

Excitations considered in the previous section constitute a part of the energy spectrum of a bosonic 
system with commensurate filling N = nL. They should not be confused with particle and hole excitations 
which arise, if we add or remove one particle from the system. The states corresponding to these particle and 
hole excitations, which are the ground states of the system with the total number of particles N = nL ± 1, 
in the lowest order of the perturbation theory are given by 




j=o 


tk 


n ± l,n,... ,n) . 

L-l 


(280) 


Up to the third order in J/U the energies of the states lip^Ld^-i) at = 0 in a d-dimensional lattice 
are |^[52] 
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4Z 


^2 


ZJ 


u 


= n-l + —n+ — (u + 1) 


U 


ZJ 


u 


+ ( ^ ) n{n+l) 


25n + ll 2n + l 
2n + l-—— + 2 


4Z 


^2 


2Z 

+ o( J") , 

5n + 1 
^ 2^ 

+ 0{J^) , 


(281) 


(282) 


where is determined by Eq. (253). 

In a one-dimensional lattice with unit filling the dependences of the energies of the particle and hole 
excitations on K up to the sixth order in J/U are given by |250| 
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Figure 24: (color online) Boundaries of the Mott-insulator phase via strong-coupling expansion up to the third order in J/U 
for d = 1 (solid lines), 2 (dashed lines), 3 (dotted lines). The gray region corresponds to = 0, and the two green regions 
correspond to the Mott insulator with n = 1,2 in one dimension. 


Eqs. (283), (284) are consistent with Eqs. (281), (282) in the same orders of J/U. 


Eqs. (281), (282), (283), (284) allow to obtain perturbative expansions for the charge gap (88). In order 


to make a comparison with the neutral gap in Eq. (279), we set d = 1 in Eqs. (281), (282) and take into 


account only the terms up to the second order in J/U. This gives 


T7 - 1 “ 2—(2?t. + 1) + ^ —^ (2r 


U 


+ 2n + 1) 


(285) 


The first-order term is the same as in Eq. (279). The second-order term is different for n > 2 and only for 
n = 1 it becomes the same. 

7 . 4 . Phase diagram 

From the energies of particle and hole excitations K — 0 one can obtain the boundaries of the regions 
in the (/i, J) plane corresponding to the commensurate fillings. They are determined by Eq. (86) and. 


therefore, are readily obtained from Eqs. (281), ( 282[ ) up to the third order of J/U for arbitrary n and d 
or from Eqs. (283), (284) up to the 6th order of J/U for n = 1 and d = 1. The region ^ < /i+(0) = —2dJ 
corresponds to V = 0. 


Eqs. (281), (282) as well as all other perturbative results are valid only in some interval of J S [0, Jc], 
where Jc is the smallest value of J at which and fJ-+{n) become equal. In this interval the difference 

/i+(n) — remains finite and positive which means that the compressibility vanishes. Jc depends on 

the order up to which the calculations are performed and provides an estimate of the transition point from 
the Mott-insulator to the superfluid. 

Perturbative results for the superfluid stiffness can be obtained from the calculations of the ground- 
state energies in the presence of the Peierls factors. In this manner we get that for N = nL‘^ the superfluid 
stiffness = 0 in all orders of J/U but the case of an incommensurate filling is different. For instance, 
in a one-dimensional lattice with N = nL ± 1, in the lowest orders of J/U we obtain 


fnL+l — 
fnL-1 = 


1 


nL + 1 

1 

nL — 1 


n -I- 1 -I- 4^n(n -|- 1) -|- ... 


n + 4:^n{n -|- 1) -|- . .. 


(286) 

(287) 


In the hard-core limit, J/U 0 and n = 1, Eq. (287) gives the same result as Eq. (236) for N = L — 1. 
Although tend to zero in the thermodynamic limit, they do not vanish in lattices of finite size. This 

is an indication that the MI lobes are surrounded by the superfluid. 


The phase diagram obtained in the third order of the strong-coupling expansion according to Eqs. (281), (282) 


is shown in Fig. 24 for different dimensions d = 1,2,3 (see also Ref. |Ml [hS]!. The size of the insulating 
regions decreases with the number of dimensions d and with the filling factor n. 
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8. Critical properties of the Bose-Hubbard model and the superfluid—Mott-insulator transition 


From the analytical results discussed in the previous sections, it becomes clear that the ground-state 
properties in the case of a commensurate filling undergo qualitative changes under variation of the ratio 
J/U. The superfiuid stiffness is equal to one for the ideal Bose gas but vanishes identically in the limit 
of strong interactions. The compressibility also vanishes for strong interactions but becomes divergent in 
the limit of vanishing interactions. Drastic changes are also observed in the excitation spectrum which is 
gapless in the case of the ideal gas but becomes gapped in the strongly interacting regime. In the case of 
an incommensurate filling, the superfiuid stiffness is always finite, the compressibility does not vanish and 
there is no gap in the excitation spectrum. 

All these characteristic features are manifestations of the SF-MI transition. Since the transition takes 
place at zero temperature, this is an example of a quantum phase transition which is driven by quantum 
fluctuations in contrast to classical phase transitions driven by thermal fluctuations. According to the 
classification adopted in the theory of critical phenomena, the transition at fixed commensurate filling 
which is controlled only by J/U and the transition controlled by the variation of the filling factor from the 
incommensurate to commensurate belong to different universality classes. 

The basic idea is that near the critical point, the spatial correlation length as well as the correlation 
time ~ where z is the dynamical critical exponent, diverge. The dependence of ^x on the distance 6 
from the critical point, which is defined as i5 = | J — Jc\/U or S = \fj, — fJ,c\/U for the two different types of 
transition, as well as the value of 2 ; do not depend on the microscopic details of the Hamiltonian and are 
entirely determined by the symmetries and system dimensionality d. The divergence of the correlation time 
leads to the vanishing energy gap 

A er' C" ■ ( 288 ) 


Dimensional analysis shows also that near the transition point the superfiuid stiffness and the compressibility 
are determined by U and U as HT] 

fs - ~ , ( 289 ) 

« ~ er" ■ (290) 


In a finite systems of L‘^ sites, the phase transition becomes a crossover. The correlation length ^x is 
limited by L and Eqs. (2881, (2891, (2901 have to be generalized. For instance, for the superfiuid stiffness, 
we have 

fs=L^-''-^<^siL/^x) , (291) 


where is a universal scaling function for this particular quantity. The power of L in front of it is the same 
as the power of ^x in the second part of Eq. (2891. The relation (2911 implies that the dependences of fs 
calculated for different system sizes L and multiplied by ];jd.+z -2 gj^ould intersect at one point which gives 
an estimation of the transition point. The expressions of the form (2911 with other powers of L also exist 
for other quantities and constitute the basis of the finite-size scaling analysis which is often used in practice 
for the calculations of the phase diagrams. 


8.1. Transition at commensurate fillings 

Now we turn to the transition at commensurate fillings which belongs to the universality class of the 
{d + l)-dimensional XY-model and characterized by z = 1 in all dimensions. It has lower critical dimension 
d = 1 and upper critical dimension d = 3. 

For d = 1, the transition is of the Berezinskii-Kosterlitz-Thouless type |275H278] . This means that the 
correlation length has an exponential divergence near the transition point 


^x ~ exp 


const 


J ^ Jq 


(292) 


and according to Eq. (288) this leads to the exponentially small energy gap. This behavior makes the study 


of the transition by computing the energy gap rather difficult and requires large system sizes. On the other 
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hand, the superfluid stiffness as well as the compressibility do not depend on [see, Eqs. (289), (290)] and, 
therefore, they do not need to be rescaled in the finite-size scaling analysis. 

The low-energy physics of the one-dimensional Bose-Hubbard model in the SF phase can be also de¬ 
scribed by the effective harmonic-fluid approach [171 127911280] and all the properties are determined by the 
Tomonaga-Luttinger parameter [171 12441128111282] 


/Ctl = (2J7r^/^(h^)K) 


- 1/2 


0 < AIxl ^ 1 


(293) 


This parameter vanishes for the ideal Bose gas and equals to one for hard-core bosons. The parameter /Ctl 
determines asymptotic behavior of the correlation functions in the superfluid phase (see sections 9.6] 9.7| ) 
and at the transition point SF-MI for the commensurate Ailing it takes the value /Ctl = = 1/2 [471 

MI phase in one dimension possesses a nonlocal string order described by the correlator [19411199] 


i 

oU() = {X{h'), 


e.'=i 


(294) 


where is the parity operator In the case J = 0, the ground state is a product of Fock states 

and Op{£) = 1 for any £. If J/U increases, the ground state contains contributions from the particle- 
hole pairs created at different distance s. As long as the positions of all particles and holes are within the 
range covered by the string correlator (294), Op{£) contains only positive contributions. However, it may 
happen, for instance, that the particle form one pair will be within the range of the correlator (294) but the 
corresponding hole not. This will give negative contribution and Op{£) will decrease. At the transition to 
the SF and above the critical point {J/U)c, the pairs are completely deconflned resulting in random positive 
and negative contributions in Fq. (294) and Op{£) vanishes. 

In the dimensions larger than one, the correlation length has a power-law dependence near the critical 
point: 

~ S-'' , (295) 


which leads to the power-law dependences of the energy gap, superfluid stiffness and compressibility. For 
d = 2, the critical exponent v « 2/3 (see, e.g., [284] and references therein). For d > 3, it takes the 
mean-field value i/ = 1/2 |47| . 

Quantum phase transitions can be viewed as fundamental changes of the ground state |4'((/)) of a many- 
body Hamiltonian H(g) = Hq + gHi under small variations Sg of the control parameter g near the critical 
point g = gc jlE]- Quantitatively, this can be described by the ground-state fidelity defined as a scalar 
product of the two ground states [285] 


$(5, 6g) = (T(5 - 6g/2)\^{g + Sg/2)) (296) 

which is expected to exhibit a sharp drop at g^- Since |dg| is small, the fidelity can be expanded as 

4>(g,Jg) = l-x(5)|+0(d4), (297) 

where xid) is a fidelity susceptibility [286] . The aforementioned drop of ^{g,6g) should be accompanied by 
a divergence of xid) which can be employed to detect the critical point [28711288] . This concept borrowed 
from the quantum information theory is appealing because it does not require any a priori identification 
of the order parameter and does not rely on the symmetries of the Hamiltonian. The scaling of fidelity $ 
and the susceptibility x near and far from the critical point is determined by the system’s dimensionality d 
and the critical exponent of the correlation length v [28711292] . Studies of fidelity for the one-dimensional 
Bose-Hubbard model of finite size are presented in Refs. It was observed that the fidelity is very 


Very often in the literature the Tomonaga-Luttinger parameter is defined as the inverse of that (Vtl 1/^tl)- 
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sensitive to the boundary conditions and its minimum is shifted from the critical point to the Mott-insulator 
side, where the changes of the ground state under variation oi g = J/U are more rapid than in the superfluid 
phase. 

The quantum phase transitions manifest themselves also through the entanglement properties of the 
system UMlIMl. This is often quantified by the von Neumann entanglement entropy 

5'a =-Tr^ (pAlnpA) , (298) 


where pA is the reduced density matrix of a subsystem A of the whole system. If the state of the whole 
system is a tensor product of the states of the subsystem A and of the remainder, Sa vanishes. In the case 
of one-dimensional Bose-Hubbard model, the SF phase is described by the conformal (i.e., relativistic and 
massless) held theory with the conformal anomaly number (central charge) c = 1. In this regime and under 
periodic boundary conditions, the von Neumann entropy of a contiguous block of La lattice sites {La < L) 
is given by 


c (L . TrLA\ , 
= - In 1 — sm I + Cl 


(299) 


where c'l is a non-universal constant. In the case of open boundary conditions, Eq. (2991 is replaced by 


c f2L . tiLa\ 


+ c'l + Sn , (300) 


where Sq, is the boundary entropy |299] . In the MI phase, the von Neumann block entropy scales as 

--ln(^a;/a) , (301) 

where the correlation length is assumed to obey inequality 1 <C ^a;/a «C La, and a is the lattice spacing. 


Eqs. (2991-(301) were used in Refs. |2501l300ll301| in order to extract the critical value of the one-dimensional 
Bose-Hubbard model from the numerical data. In Ref. [MU the critical point in one dimension was obtained 
using the concept of bipartite fluctuations |8()8Hfl0fl] through the analysis of the particle number fluctuations 
in the subsystem A: 

Ta = (^1) - {^A? , TVa = 5] ni . (302) 

IgA 

This quantity has similar scaling properties as the von Neumann entropy and allows quite efficient calcula¬ 
tions of the Tomonaga-Luttinger parameter A^tl- 

Exact calculation of the critical values of J/U in different dimensions and for different fillings is a 
challenging problem and various methods have been applied in order to achieve this goal. A brief overview of 
these extensive studies in one dimension at unit Ailing was given in Ref. (MU. These results are summarized 
in Table [2 which includes also some other missing as well as more recent references. Although the numerical 
data are quite different, all methods give much larger values than the prediction of the mean-field theory 
{J/U)c ~ 0.086 [see Eq. (332l]. Since the results obtained by exact diagonalization suffer only from the 
finite-size effects, they can be considered as a lower estimate of {J/U)c- Taking into consideration also 
the most recent studies for larger systems, one can accept that {J/U)c ~ 0.3. Another observation is that 
analytical methods have a tendency to underestimate the critical value {J/U)c- 

The critical values of J/U in one dimension and in the case of unit Ailing were measured first in experi¬ 
ments with Cs atoms using the lattice modulation spectroscopy |B]. For the lattice depths Vq/Ar = 6 ... 10 
corresponding to the tight-binding regime, the experimental values of {J/U)c appeared to be less than 0.26 
calculated in Ref. |165| . This systematic underestimation was attributed to the presence of the harmonic 
trap that leads to the spatial inhomogeneity and finite-size effects. Similar experiment was performed very 
recently with |323| . The transition points determined from the measurements of the critical momentum 
for the occurrence of a dynamical instability are in good agreement with {J/U)c ~ 0.297. 

For larger integer fillings, exact diagonalization cannot be applied anymore due to very limited systems 
size but other methods still can be used and some of the references listed in Table [^reported also estimations 
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iJ/U)c 

method 

quantity 

L 

year 

Ref. 

0.215 ±0.01 

QMC 

energy gap 

up to 32 

1990 

12QZI 

l/(2v^) « 0.2887 

BA 

superfluid stiffness 

oo 

1991 


0.215 

RSRG 

flxed point 

oo 

1992 

|HnR| 

0.215 

SCE 

energy gap 

oo 

1994 

m 

0.275 ±0.005 

ED 

energy gap 

up to 9 

1994 

IM] 

0.22 ±0.02 

ED 

energy gap 

up to 11 

1995 


0.298 ±0.002 

DMRG 

energy gap 

up to 70 

1996 

|310| 

0.304 ±0.002 

ED±RG 

superfluid stiffness 

up to 12 

1996 

m\ 

0.300 ±0.005 

QMC 

energy gap 

up to 50 

1996 

I3TT1 

0.265 

SCE±RG 

energy gap 

OO 

1996 

isa 

0.25 

TDVP 

energy gap 

oo 

1998 

pra 

0.277 ±0.01 

DMRG 

Fa 

up to 76 

1998 

|313| 

0.26 ±0.01 

SCE±PA 

energy gap 

OO 

1999 

|251| 

0.260 ±0.005 

DMRG 

superfluid stiffness 

up to 50 

1999 

11 ^ 

0.297 ±0.01 

DMRG 

Fa 

up to 1024 

2000 

m 

0.283 ±0.005 

ED 

superfluid stiffness 

up to 14 

2004 

uni 

0.305 ± 0.004 

QMC 

Fa 

128 

2005 

lasi 

0.257 ±0.001 

ED 

ground-state fldelity 

up to 12 

2007 


0.238 ±0.011 

GFMC 

structure factor 

up to 150 

2007 

|316| 

0.204 ±0.004 

VMC 

structure factor 

up to 150 

2007 

m 

0.303 ±0.009 

DMRG 

energy gap 

up to 80 

2008 

uni 

0.2975 ± 0.0005 

TEBD 

Fa 

OO 

2008 

m 

0.29... 0.30 

DMRG 

von Neumann entropy 

up to 1024 

2008 

l^nni 

0.305 ±0.001 

DMRG 

Fn 

up to 1024 

2011 

|M9j 

0.319 ±0.001 

TEBD 

energy splitting 

up to 48 

2011 

|H5ni 

0.295... 0.320 

DMRG 

string correlator 

216 

2011 

IMj 

0.2989 ± 0.0002 

DMRG 

bipartite fluctuations 

up to 256 

2012 

pna 

0.2885 ± 0.0001 

NBA 

superfluid stiffness 

up to 1400 

2012 

|321] 

0.30 ±0.01 

TEBD 

von Neumann entropy 

OO 

2012 

I3UT1 

0.305 ±0.003 

DMRG 

von Neumann entropy 

up to 64 

2012 

12501 

0.3050 ±0.0001 

DMRG 

energy gap 

up to 700 

2013 

pM] 

0.270 ±0.008 

TEBD 

ground-state fldelity 

up to 64 

2014 

[295] 

0.289 ±0.008 

TEBD 

fldelity susceptibility 

up to 64 

2014 

[295] 

0.286 ±0.005 

ED 

energy gap 

up to 12 

2015 

j322j 


Table 1: The critical value of J/[/ for the MI-SF phase transition in one dimension at unit filling obtained by the analysis of 
different physical quantities in the lattices of size L employing different method such as quantum Monte Carlo (QMC), Bethe 
Ansatz (BA), numerical solution of Bethe equations (NBA), real-space renormalization group (RSRG), strong-coupling expan¬ 
sion (SCE), exact diagonalization (ED), density-matrix renormalization group (DMRG), Fade analysis (PA), time-evolving 
block decimation (TEBD), mean-field theory based on the time-dependent variational principle (TDVP), Green’s function 
Monte Carlo (GFMC), variational Monte Carlo (VMC). 
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(J/U)c 

method 

quantity 

L 

year 

Ref. 

0.061 ±0.003 

PIMG 

superfluid stiffness 

up to 8 

1991 

13251 

0.0564 

RSRG 

fixed point 

oo 

1992 

|HM| 

0.0585 

SGE±RG 

energy gap 

oo 

1996 

m 

0.0625 

TDVP 

energy gap 

oo 

1998 


0.05974 ±0.00004 

SGE±PA 

energy gap 

oo 

1999 

|2BT| 

0.05963 ±0.00001 

SSE 

superfluid stiffness 

up to 20 

2005 

1326] 

0.0485 ± 0.0005 

VMG 

structure factor 

up to 30 

2008 

1327] 

0.0588 ± 0.0007 

GFMG 

structure factor 

16 

2008 

IB27| 

0.05974 ±0.00003 

WA 

energy gap 

up to 80 

2008 

m\ 

0.05909 

MEP±SGE 

susceptibility 

OO 

2009 

ca 

0.067 

VGA 

energy gap 

oo 

2010 

132^ 

0.060 

NPRG 

superfluid stiffness 

oo 

2011 

|330l 

0.055 

POA 

superfluid stiffness 

oo 

2011 

PH 


Table 2: The critical value of J/C/ for the MI-SF phase transition in two dimensions at unit filling obtained by the anal¬ 
ysis of different physical quantities in the lattices of linear size L employing different method such as path integral Monte 
Carlo (PIMC), stochastic series expansion (SSE), worm algorithm (WA), variational Monte Carlo (VMC), Green’s func¬ 
tion Monte Carlo (GFMC), real-space renormalization group (RSRG), strong-coupling expansion (SCE), Fade analysis (PA), 
mean-field theory based on the time-dependent variational principle (TDVP), nonperturbative renormalization group (NPRG), 
method of effective potential (MEP), projection-operator approach (POA). 


of {J/U)c in one dimension for (fi^) = 2,3. Analytical studies within the SCE up to the third order in 
combination with the scaling theory found (J/t/)c = 0.155 for {rif} = 2 and {J/U)c = 0.111 for {hi) = 
3 [S5]. Slightly lower value of {J/U)c — 0.008 for {fif} = 3 was obtained in the mean-field theory based 
on the TDVP |312l 1324] . Numerical study of the correlation function Fa{s) with the aid of the TEBD 
method designed for infinite systems gave {J/U)c = 0.175 ± 0.002 for {n() = 2 |318| which is very close 
to {J/U)c = 0.180 ± 0.001 obtained from the DMRG calculations of the correlation function F„(s) in the 
systems up to L = 128 |319| as well as to ( J/t/)c = 0.179 ± 0.007 resulting from the DMRG calculations of 
the von Neumann entropy in the systems up to L = 64 pso] . Recent DMRG calculations of the energy gap 
for the systems up to L = 250 lead to (J/{7)c = 0.1790±0.0003 for (fii) = 2 and {J/U)c = 0.12697±0.00003 
for {nt) = 3 |294| . 

In Ref. |320| . the critical points in one dimension were calculated for the fillings up to {fit) = 1000 from 
the energy splitting caused by the tunneling between two states with macroscopically distinct currents using 
the TEBD method. The numerical data were well fitted by the function 


= d{ni) (a + b{ne) °) , 


(303) 


with the coefficients a = 2.16, b — 0.97, c = 2.13. 

In higher dimensions, the arsenal of the methods is restricted because exact diagonalization and the 
DMRG cannot be applied. Nevertheless, QMG works very well and approximate analytical methods can be 
also used. The results for two-dimensional square lattices at unit filling are summarized in Table It is 
interesting to note that one (semi)analytical and two different numerical methods used in Refs. |2511l326ll328| 
give essentially the same value of (J/C/)c ~ 0.0597. The method of the effective potential, developed in 
Refs. |79[ IBP] as a combination of the mean-field theory and high-order SGE, as well as the nonperturbative 
renormalization group approach |330| lead to almost the same results which are considerably larger than 
{J/U)c ~ 0.0429 predicted by the standard mean-field theory [see Eq. (332l]. 


For larger fillings in two dimensions, to the best of our knowledge, no exact numerical results were 
published but some approximate analytical results are available. The third-order SGE after extrapolation 
to the infinite order gives {J/U)c = 0.0345 for (hi) = 2 and {J/U)c — 0.0245 for (hi) = 3 [S5], while the 
VGA results in {J/U)c = 0.038 for (hi) = 2 |329| . 

In three-dimensional cubic lattices, the critical value of J/U was calculated by the extrapolation of the 
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third-order SCE to the infinite order which gives {J/U)c = 0.0337 ± 0.0027, 0.0200 ± 0.0013, and 0.0140 ± 
0.0007 for (hi) = 1, 2, and 3, respectively [52]. The numerical value for (hi) = 1 is in excellent agreement 
with QMC calculations of the energy gap in the lattices with the linear sizes up to L = 20, which give 
{J/U)c = 0.03408±0.00002 |204| . The NPRG approach yields almost the same result {J/U)c = 0.0339 [330| 
and the POA developed in Refs. |331L 1332] allows to reproduce the results of the QMC calculations with 
the accuracy ~ 0.05%. Although the mean-field theory is expected to work better in higher dimensions, 
Eq. (332) provides an estimate (J/C/)c ~ 0.0286 which is again noticeably lower than the exact numerical 
values. VMC calculations of the static structure factor for systems with linear sizes up to L = 12 gave 
{J/U)c ~ 0.0278 |316| . Recently developed self-consistent standard basis operator approach, which goes 
significantly beyond the mean-field theory, yields {J/U)c ~ 0.03356 for (hi) = 1 and {J/U)c ~ 0.0185 for 
(hi) = 2 |3Mj. 

Theoretical predictions for the critical value { J/U)c in three dimensions were tested in experiments with 
Cs atoms analyzing the quasi-momentum distributions in the time-of-flight images pi] . It was observed that 
the width of the central peak as a function of the amplitude of the periodic potential Vq at fixed scattering 
length Os has a kink at some value V) that is interpreted as the transition point. The experimental data 
obtained in deep lattices (Vb = 8... 18 Ar) in the case of unit filling for different values of Ug controlled 
with the aid of Feshbach resonances are in good agreement with the critical values obtained by QMC |204| 
as well as in the framework of the mean-field theory within the experimental uncertainty. 

In Refs. |79l [8fT| . the method of the effective potential combined with the high-order SCE was used for 
the calculations of {J/U)c for fillings up to (hi) = 10000 in two and three dimensions. The computed critical 
values were fitted by 

J\ ^ J \ 0 13 

^ \ f 0.13 


u 


u 


>/(«i) {{n\) + l)d'^ ' 


where (J/C/)“^ is given by Eq. (332), with the accuracy of about 1% [SD], and by 


MF 


1 


0.35 


0.39 


0.84 

~dF 


(305) 


with the accuracy of 0.15% |335| . It was pointed out |320| that the fit (303) with the coefficients (a, 6, c) = 
(5.80, 2.66, 2.19) and (a, 6, c) = (6.70, 3.08, 2.18) in two and three dimensions, respectively, gives also a good 
approximation for the transition points. 


8.2. Generic transition 

The transition governed by the variation of the filling factor (or chemical potential) belongs to the mean- 
field universality class (Caussian model). The behavior of the spatial correlation length near the critical 
point is described by Eq. (295) and the critical exponents are z = 2 and v = 1/2 in all dimensions This 
transition has the upper critical dimension d = 2. In one dimension, the Tomonaga-Luttinger parameter 
takes its universal value /Ctl = tl = 1 at the transition points |24411281] . 

The critical points of the generic transition are described by two continuous lines /i±(J) parametrized 
by the integer filling factors (hi) which cross at the critical point of the commensurate transition and form 
the boundaries of the insulating regions in the plane (MI lobes). Within the third-order SCE, the 

boundaries /i±(J) are given by Eqs. (86), (281), ( |282[ ) and the corresponding phase diagrams are shown 
In Figs. I 


in Fig. In Figs. [^ we show the exact phase diagrams worked out by numerical and semi- 

analytical methods. The comparison with Fig.[^shows that few lowest orders of strong-coupling expansion 
appear to be sufficient in order to reproduce the topology of the phase diagram in three and two dimensions 


(Figs. 25 26) except the tips of the lobes. Higher-order calculations [2511 1252] and extrapolation to the 
infinite order [52 1 IM] I25n 1252] allow to reach perfect agreement with exact numerical results. QMC data 
for the filling factor (hi) = 1 in two |328| and three |204| dimensions were also reproduced with the high 
accuracy using the MEP |79[ I5IJII555] . the B-DMFT |851[55] . and the NPRC approach |330| . The POA gives 
also a high accuracy in three dimensions 155111552] . 
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Figure 25: (Color online). Ground-state phase diagram of the Bose-Hubbard model in three dimensions showing MI lobes 
with (hi) = 1,2 surrounded by the SF. (Adapted with permission from Ref. [86]). 



Figure 26: (Color online) Ground-state phase diagram of the Bose-Hubbard model in two dimensions showing MI lobes with 
(hi) = 1, 2 surrounded by the SF. (Adapted with permission from Ref. ESI). 


In one dimension, the situation is more complicated due to the fact that the critical behavior of the 
system is of the Berezinskii-Kosterlitz-Thouless type. The shape of the boundaries fj,± separating the MI 
from the SF is qualitatively different (see Fig. 271. Starting from certain values of J/U, the lowest boundary 
fj,- bends down which implies that in some interval of /i the MI phase is reentrant, i.e., increasing the 
tunneling parameter J one returns to the MI. This feature becomes visible in the SCE calculations up to 
the 12th order |251L 1313] , 


8.3. Finite temperature 

At finite temperature, thermal fluctuations give rise to the normal phase which appears on the phase 
diagram in addition to the SF and MI phases. Fig. shows the finite-temperature phase diagram of a 
two-dimensional system worked out by QMC simulations [557] . and in three dimensions the topology should 
remain the same. Compared to T = 0, the boundary separating SF from the insulating phases is shifted 
towards larger values of J/U. Due to the fact that the compressibility never vanishes at finite temperature, 
there is no drastic difference between the MI and normal gas. On the other hand, if the compressibility 
K is small enough, the system can be still considered as a MI. In the example for two-dimensional system 
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Figure 27: Ground-state phase diagram of the Bose-Hubbard model in one dimension showing MI lobes with (h^) = 1,2 
calculated by the DMRG method for L < 128. (Adapted with permission from Ref. EISl). 


shown in Fig. |28| the crossover line between the MI and normal gas was determined from the requirement 
kU < 0.04 [337|. 





(b) 


15 li/j 


Figure 28: (color online), (a) Finite-temperature phase diagram for the homogeneous Bose-Hubbard model in two dimensions 
in the {fi/U,t/U) plane for k^TjV = 0.02, 0.04 and 0.08. The lines that demarcate the SF and N is a phase boundary, and 
the lines that demarcate MI and N is a crossover. At finite temperature, normal phase regions appear between MI and SF. 
These normal regions are bigger for higher temperature. (Adapted with permission from Ref. I337l . ©2011 American Physical 
Society.) (b) Critical temperature in two dimensions at filling {hi) = 1. Circles are simulation results and the dotted line is to 
guide the eye. Dashed lines are analytical results for the weakly interacting gas | |308[ l and for the strongly interacting gas near 
the quantum critical point ( |307[ |. (Adapted with permission from Ref. 1^ . ©2008 American Physical Society.) 

On the other hand,, the compressibility in the insulating phase has a “thermally activated" form, i.e., 
k{T) ~ exp(—A/^B^r) with a finite energy gap A (see, e.g.. Ref. |338| 1. This suggests an alternative 
specification of the crossover line between the normal gas and the MI from the condition A = A:bT |338F 

EH]. 

Figs. mb) , show the dependences of the critical temperature T^. of the superfluid-normal transition 
in two and three dimensions for an integer filling (hi) = 1 |20411328] . Near the quantum critical point of the 
SF-MI transition (U/J)c at zero temperature, the critical temperature is related to the superfluid stiffness 

as mmsiEia 
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Figure 29: Critical temperature in three dimensions at filling {hi) = 1. Circles are QMC simulation results and the line is to 
guide the eye. / J = 5.591 is the critical temperature of the ideal Bose gas with the tight binding dispersion relation. 

(Adapted with permission from Ref. 1^ . ©2007 American Physical Society.) 


Taking into account Eqs. (2891, (2951, we obtain 


k^Tc 

J 



(307) 


where z = 1, and i/ is the critical exponent of the correlation length of the (d + l)-dimensional XY-model. 
In two dimensions, QMC calculations gave A = 0.49 ± 0.02 |328| . 

In the weakly interacting regime, the dependence of Tc in two and three dimensions is qualitatively 
different. In two dimensions, the transition is of the Berezinskii-Kosterlitz-Thouless type and for small U/J 
the critical temperature is given by [32811343| 


J 


47r(fii) 
ln(2^J/C/) ’ 


C = 380 ± 3 . 


(308) 


This equation shows that in the limit of vanishing interaction, T™ vanishes. In three dimensions, the 
critical temperature in the limit of vanishing interaction tends to a finite value (see section 6.1.31. The 
critical temperature of the superfluid-insulator transition experimentally measured in three-dimensional 
optical lattice 


shows satisfactory agreement with the QMC data presented in Fig. 29 


Quantum critical behavior was observed in experiments with Cs atoms in a 2D optical lattice at finite 
temperature near the normal-to-superfluid transition m- In the zero-temperature limit, it connects to 
the vacuum-to-superfluid transition, where vacuum can be viewed as a MI with zero occupation number. 
On the basis of in situ density measurements, the equation of state p{^, T) of the sample was determined, 
which gave the values of the critical exponents z = 2 .2'^^'^ and v = 0.52^° °® in agreement with theoretical 
predictions. 


8.4- Criticality in confined systems 

In order to create a MI in a homogeneous lattice, it is necessary to have and absolute control over the total 
number of particles N which has to be a multiple integer of the total number of lattice sites L'^. This is hardly 
achievable with ultracold atoms. However, real experiments are performed in harmonic traps which allow 
to create MI in restricted spatial regions. The presence of harmonic confinement leads to inhomogeneous 
spatial distribution of (hi) in the ground state |147l 114811344H346| . see Fig. 30 For sufficiently small total 
number of particles, the density profile is smooth having the form of the inverted confining potential. With 
the increase of the particle number a plateau with a local filling of one boson per site develops in the central 
part of the trap, provided that the number of particles exceeds some critical value which depends on the 
system parameters, similarly to the case of hard-core bosons considered in section |6.4.3| In the example 
shown in Fig.j^this critical number is about 30. Within the plateau the local compressibility is small which 
is considered as an indicator of the MI region. Further increase of the particle number leads first to the 
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Figure 30: (a) Spatial distribution of the mean occupation numbers (h^) in a one-dimensional lattice in the presence of harmonic 

confinement for increasing total number of bosons N calculated by QMC. The parameters are V^jJ = 0.008, U/ J = 8, L = 100. 
At low N, the system is in a SF phase. MI plateaus appear as N is increased, (b) Profiles of the local compressibility ki 
(solid lines) associated with corresponding mean particle-number distributions (circles) for the same parameters as in (a) and 
for fixed TV. is very small when {hi) = 1. (Adapted with permission from Ref. Ez]. ©2002 American Physical Society.) 


broadening of the plateau but then to the formation of a compressible region with local fillings larger than 
one in the center of the trap. Increasing the particle number further, one can produce a second MI region 
leading to the so-called wedding-cake (shell) structure of the density distribution, which can be interpreted 
as a coexistence of the SF and MI phases. 

In the presence of confining potentials it becomes more reasonable to consider state diagrams instead of 
phase diagrams. It turns out that the state diagrams are determined by two parameters, the characteristic 
density p and J/U, and almost independent on the values of V^j J |:i45l KUb] . In the case of harmonic trap, 
p = Na‘^{VT/J)‘^^'^ |206l 13451 . It was found that the largest values of J/U that support the local insulator 
with {flu) = 1 are 0.18 in one dimension and 0.0575 in two dimensions |346| . The latter is very close to the 
critical value of the SF-MI transition in a homogeneous 2D system, but the former is substantially lower 
than the corresponding results in ID (see Tables mi)- 

Since diverging length scales cannot appear in confined systems, it was debated whether quantum criti¬ 
cality can be observed in experiments pRlITTmiMillEK] . In earlier papers, the experimental observation 
of the MI phase was interpreted as a crossover |148| . Later it was demonstrated that the critical behavior 
of trapped systems can be cast in the form of trap-size scaling in analogy to the finite-size scaling theory 
for homogeneous systems, which shows that at criticality the spatial scale depends on the trap size R as 
~ R^ with an additional trap critical exponent e fmiiM7| . The value of 6 can be obtained using the 
renormalization-group analysis. In the case of power-law potentials with the exponent p {p = 2 corresponds 
to the harmonic trap) and for Bose-Hubbard model. 


p+l/u 


(309) 


where v is the critical exponent of the correlation length of the homogeneous system m- It was shown 
that violation of the local-density approximation is a good indication of the critical behavior, which can 
be observed from the measurements of in situ density profiles and quasi-momentum distribution |348| . 
Experiments in two-dimensional lattices confirm this prediction and the critical values of J/U extracted from 
the measurements of the fraction of particles with zero momentum are consistent with QMC calculations 
for trapped systems |23| . 


9. Exact numerical results 

In this section we present exact results for the Bose-Hubbard model in different dimensions d= 1,2,3 for 
finite systems across the quantum critical point. They are obtained by different numerical methods such as 
exact diagonalization, density-matrix renormalization group (DMRG) and quantum Monte Carlo (QMC). 
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Exact diagonalization is performed in the complete Hilbert space of the Bose-Hubbard Hamiltonian (541 


and allows to calculate in principle all the eigenvalues and all the eigenstates. However, due to the rapid 


growth of the Hilbert space with the number of particles and the number of lattice sites [see Eq. (611] the 
calculations are limited to the systems of small sizes. In homogeneous lattices under periodic boundary 
conditions the dimension of the Hamiltonian matrices can be reduced approximately by a factor L‘^ using 


the basis of the eigenstates of the momentum operator that are given by Eq. (641. However, this does not 


allow to increase considerably the size of the systems. On the other hand, typical eigenstates of quantum 
many-body systems occupy only a small part of the Hilbert space. This idea is used in the DMRG method 
which is very efficient in one dimension and allows to treat large systems. In higher dimensions, QMC 
methods based on the stochastic sampling of the complete Hilbert space appear to be superior. 

All these methods as well as some others are implemented in the ALPS package which is freely avail¬ 
able [Ml- It is not difficult to install it and it was already used in many publications (see a list of references 
on the homepage of the ALPS project). The results of the DMRG and QMG calculations reported in this 
section are obtained with the aid of the package. However, comparison of our own implementation of the 
exact diagonalization with that of ALPS showed that the latter is inefficient and requires too much computer 
memory, therefore, we used our own exact diagonalization code. Before we turn to the discussion of results 
we would like to make some notes about exact diagonalization. 

9.1. Remarks on exact diagonalization 

In simplest situations like in the case of two atoms considered in section [63] one can write down explicitly 
all the basis states and all the matrix elements of the Hamiltonian. However, for arbitrary number of 
atoms and lattice sites this is not possible and we need efficient computer algorithms to deal with this. In 
order to generate the Hamiltonian matrix, we have to generate sequentially all the basis states, act by the 
Hamiltonian on each of the states and determine the numbers of the resulting states. This leads us in the 


field of combinatorics because the basis states (61) are nothing but compositions of an integer N into L non¬ 


negative parts. There are efficient algorithms that fulfill the following tasks for given N and L: (1) generate 
sequentially all compositions; (2) determine the number of a given composition; (3) generate a composition 
with a given number. The tasks (2) and (3) are called ranking and unranking, respectively. All these tasks 
are implemented in the combinatorial package SELEGT [35011351] and can be use for matrix generation as 
well as for calculations of the observables after the diagonalization. However, the package does not take 
into account the spatial symmetries like translational invariance or discrete reflections and some additional 
programming is required in order to use those. Alternatively, one can use hashing techniques |69l 1352] . 

If we want to make a profit from the translational invariance, we have to use the basis states (65l. The 


interaction part of the Bose-Hubbard Hamiltonian creates only diagonal matrix elements which is simple. 
The hopping part of the Hamiltonian is non-diagonal and in order to understand how to create corresponding 
matrix it is sufficient to consider the action of an operator which describes hopping 

“from right to the left", on a basis state ([6^. We note that 


aja^_^jnr) = W(nr,£ + 1) nr,£+i |nr,) 


(310) 


where rur^ is an integer which is uniquely determined by T and i. Then for fixed K it is easy to show that 


L _ . - 

^ \/ {nr,£ + 1) nr.£+iA [nicr^) • 

7^1 V 


(311) 


One has to keep in mind that the states [njcr^) in Eq. (311 \ are not necessarily distinct and may coincide 


with Injcr)- Galculating all the matrix elements of the operator one can construct the full Hamiltonian 
matrix. 


9.2. (p-, J) diagram 

We consider first the boundaries fJ,±{N) of the regions in the (/i, J) plane corresponding to different 
total particle numbers N. They are determined by Eq. (86l. The results of numerical calculations for one¬ 
dimensional lattices obtained by exact diagonalization are shown in Fig. |3H and the qualitative behavior 
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Figure 31: (a) Regions with the occupation numbers N = 0,1,..., 15 (from bottom to the top) in a one-dimensional lattice 

of L = 10 sites, (b) The red lines are the boundaries of the region with N = L calculated for L = 2, 3,..., 13 by exact 
diagonalization. The black lines are the boundaries from the DMRG calculations for 128 sites. See also Refs. |293l 131^ I353| . 
(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 




Figure 32: Superfluid stiffness obtained by numerical diagonalization of the Bose-Hubbard Hamiltonian for L = 10 in the 
presence of Peierls factors with rj} = 0.01. (a) f‘^ as a function of N for J/U = 0 (red), 0.08 (green), 0.16 (blue), 0.24 (magenta), 
0.32 (brown), 0.5 (black), (b) as a function of J/U for N ranging from 1 (/^ = 1) to 10 from top to the bottom. Dashed 
line is the result of the second-order perturbation theory in J/U [Eq. | |287| ]. (For interpretation of the references to colour in 
this figure legend, the reader is referred to the web version of this article.) 


remains the same in higher dimensions. Due to the finite system size, the regions of different occupation 
numbers are always finite but their boundaries come closer to each other with the increase of the number 
of lattice sites. In the thermodynamic limit, they should densely cover the whole (/i, J) plane except some 
finite regions corresponding to the MI phase which exist only for commensurate fillings, provided that the 
ratio J/U is small enough. Due to the fact that the number of particles does not depend on /i within these 
finite regions, the compressibility (821 vanishes, and we should obtain exactly the same phase diagram as in 
Fig.l^ 


9.3. Superfluid stijjness 

The superfluid stiffness calculated according to the definition ( [M| ) in a finite one-dimensional lattice 
for different particle numbers N and different values of J/U is shown in Fig. 32 If the filling is not 
commensurate, the superfluid stiffness remains finite for all values of J/U. However, for commensurate 
fillings, finite intervals of J/U close to zero appear, where with a good numerical accuracy can be 
considered as vanishing, i.e., the system becomes an insulator. Numerical values of in the limit J/U —?► 0 


are perfectly described by Eqs. (2361, (2521 derived for hard-core bosons in one dimension. 


In the hydrodynamic regime, the superfluid stiffness is determined by Eq. (95). Together with Eq. (144) 
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Figure 33: (color online) Comparison of the exact results for the superfiuid stiffness (black solid lines, the same as in Fig. |32[ 3) 
with those obtained from Eq. | |312[ | within the framework of the hydrodynamic approach (red dashed lines). 


for the sound velocity at T = 0 derived from the sum rules for the dynamic structure factor, it gives 


/^ = 


(«i) 


(312) 


For the Tonks-Girardeau lattice gas, this expression leads to Eq. (2361. However, in general the validity of 


Eq. (312) is limited. This is demonstrated in Fig. 33 where one can see that Eq. (312) works well only for 


low fillings. For fillings close to one, there are strong deviations from the exact results at small J/U. 

9 . 4 . Energy spectrum 

As it was discussed in Sec. |3.3| the eigenstates of the homogeneous system under periodic boundary 
conditions are characterized by two indices, K and fl, where the former denotes the total momentum of N 
interacting particles in a lattice. In the case of commensurate filling, the number of eigenstates for K = 0 
is always larger than for other values of K. For the discussion of the energy spectrum, it appears to be 
convenient to start the labeling for AT = 0 with fl = 0 and for other K’s with = 1. 

The full energy spectrum calculated by exact diagonalization for = L = 11 and J/U = 0.05, which 
corresponds in the thermodynamic limit to the Mott-insulator state, is shown in Fig. [Mfa). The lowest dot 
at AT = 0 is the ground-state energy see also Fig. 


34 


b). The energies = 1,..., L — 1, form 


the lowest excitation band shown by the solid lines in Figs. 34 ’b,c). For such a small value of J/U, it does 
not overlap with the higher excitation bands and is well described by Eq. (278). An increase of the system 


size leads to more dense distribution of the points, however the structure of the lower part of the spectrum 
remains the same [compare Figs. mb) and [M|[c)]. At small momenta K, the lowest excitation branch can 
be approximated by a pseudo-relativistic form 


{EP-El°y = iA£y+v%K^ 


(313) 


with the energy gap AS and the effective velocity Ueff. 

If we increase the ratio J/U, the excitation bands start to overlap and the energy spectrum becomes 
qualitatively different (see Fig. m^) and Ref. m)- The lowest excited state is not at Kq = 0 anymore 
but at K±i = ±27r/(La) and, therefore, it is degenerate. In the thermodynamic limit, this should give 
a sound mode with the linear dispersion relation for small K. In the limit of infinite interaction and for 


commensurate fillings, the corresponding sound velocity vanishes according to Eq. (233) but remains finite 


for finite values oi J/U above the critical point |355| . In higher-dimensional lattices, the structure of the 
energy spectrum is expected to be the same, although no exact results for sufficiently large systems have 
been reported so far. 

Energy-spectrum statistics of the one-dimensional Bose-Hubbard model was studied in Refs. |3561I359] . 
While in the special cases J = 0, C/ = 0, as well as in the hard-core limit, the model is obviously integrable. 
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Figure 34: (color online) (a) Full energy spectrum and (b) its lowest part for N = L = 11, J/U = 0.05. (c) Dots: lowest part of 
the energy spectrum for N = L = 14, J/U = 0.05. Shaded region is the lowest excitation band obtained in the thermodynamic 
limit in the second order of J/U [Eq. ( |278[ l]. Solid lines show the boundaries of the lowest excitation band which is well resolved 
for this value of J/U. (d) Lowest part of the spectrum for N = L = 14, J/U = 0.2. Solid line connects the energies of the 
lowest excited states. 


at intermediate values of the parameters the integrability is lost. The distribution of spacings between the 
adjacent energy levels follow the universal Wigner-Dyson law characteristic for quantum chaotic systems 
that belong to the Gaussian orthogonal ensemble. In particular, the probabilities of small level spacings are 
strongly suppressed which is an indication of level repulsion due to the avoided crossings. 

In Figs. |35| we compare exact numerical data for the ground-state energy in different dimensions with 
the corresponding results of the strong-coupling perturbation theory (2531, (2541. We observe that the 4th 
order SCE gives already a satisfactory description up to the critical point of the SF-MI phase transition 
in all dimensions. We have also compared the results of exact diagonalization for 14 sites in ID with 
the DMRG-calculations for larger systems and did not find any noticeable discrepancy which shows that 
finite-size effects (at least for local quantities) are negligible for systems of this size. 

The comparison of the ground-state energy in two dimensions calculated by exact diagonalization for 
lattices of 3 x 3 sites with QMG data for larger lattices reveals a contribution of the finite-size effects Note 
that the numerical values for a lattice of 10 x 10 sites are better described by the SGE in the thermodynamic 
limit than those obtained for 3 x 3 sites. We also compared the results of exact diagonalization and QMG 
for 3x3 sites and found perfect agreement not only for the ground-state energy but also for other quantities 
discussed below. 

In three dimensions, exact diagonalization does not make much sense due to strong finite-size effects 
and we have to rely on QMG. The numerical data obtained for the lattice of 5 x 5 x 5 sites are in a good 
agreement with the SGE near and below the critical point. 

The energies of neutral and charge excitations in a one-dimensional lattice in the case of unit filling 
are plotted in Fig. As it was discussed in section |7.3| An and Ac coincide at least in the second-order 
of the strong-coupling expansion. Numerical results presented in Fig. confirm this prediction at small 
values of J/U but show that An and Ac become different for larger JjU . In the thermodynamic limit, both 
quantities are expected to vanish above the critical point { J/U)c- However, in a finite system, Ac remains 
almost constant for increasing J/U and A^ even grows. In addition, the latter has a point of nonanalyticity 
marked by a vertical dotted line in Fig. 36 For J/U below the point of nonanaliticity the lowest excited 
state of the system is at Kq = 0 and above this point it is at AT = K±i, i.e., we jump from one excitation 
branch to the other [see Fig.[M|(b,d)]. The growth of An with J is a finite-size effect which can be understood 
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Figure 35: Ground-state energy for unit filling in one (a), two (b), and three (c) dimensions. Solid lines are the results of 
exact diagonalization for 14 sites in ID (a) and for 3x3 sites in 2D (b). Dots are QMC data for 10 X 10 sites in 2D (b) and 
5x5x5 sites in 3D (c). Dashed lines in all panels are the results of the strong-coupling expansion up to the fourth order 
[Eq. l|253^]. Dotted line in (a) is obtained by the strong-coupling expansion up to the 14th order [Eq. ||254[||. 


looking at the limit L —oo in Eq. (1541 for the ideal Bose gas. 


9.5. Particle-number distribution 

In the present section, we discuss the probabilities p{n\ = n) to have n atoms at a lattice site in the case 
of unit filling, N = L'^. Due to translational invariance, p{n\ = n) does not depend on the site index. At 
zero temperature and in the limit J = 0, we have p{n\ = n) = Sn.i, and in the opposite limit, U = 0, the 
probabilities are given by the binomial distribution ( [T58| . For finite J and U, exact numerical results are 
shown in Fig. One can clearly see that the probabilities to have three particles or more at one lattice 
site are very small for any J/U in the case of unit filling. For small values of J/U, p{n\ = 0) and p{n\ = 2) 
are almost equal to each other, no matter what the system dimensionality is, which is a manifestation of the 
particle-hole symmetry. The exact numerical data are in good agreement with the SCE. Below {J/U)c, the 
probabilities of the occupation numbers different from one become smaller if the dimensionality is increased 
(see Fig. [37| . In Fig. [37|(b) one can see that the finite-size effects lead to the larger values of p{ni — 0) 
and p{n\ = 2) below the critical point and to the lower values of those above (J/t/)c compared to the 
thermodynamic limit. 

The particle-number fluctuations ant — \/{’^D ~ (standard deviation) at finite temperature cal¬ 

culated by exact diagonalization in one dimension are shown in Fig. |37[ a’). Comparison with the zero- 
temperature result for the same system size indicates that temperature has stronger influence at smaller 
values of J/U. 


9.6. One-body density matrix 

Now we turn to the discussion of the two-point correlation functions. In a homogeneous lattice under 
periodic boundary conditions they depend on F —li . F irst, we consider the one-body density matrix with the 


entries Fa (li, h) = As one can see in Fig. 38 it monotonically decreases with s and increases with 

J/U approaching the asymptotic limit of the idealgas (1611 for J/U —>■ oo, and the results of numerical 
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Figure 36: Energies of neutral (red lines) and charged (green line) excitations in a one dimensional lattice for unit filling. 
Solid lines - exact diagonalization for L = 13, red dashed line - strong-coupling expansion of the 6th order [Eqs. ( |283[ l, ( |284[ | 
with K = 0]. Vertical dotted line marks the point, where = EjJ. (For interpretation of the references to colour in this 
figure legend, the reader is referred to the web version of this article.) 


calculations for small J/U are in good quantitative agreement with the SCE. Higher-order SCE in one 
dimension allows to extend the region of validity of the perturbative calculations. 

In the MI phase, the presence of the energy gap in the excitation spectrum leads to the exponential 
asymptotics of the two-point correlation functions in the ground state at large distances IMU. This general 
statement is valid in any finite dimension and the corresponding length scale is inversely proportional to the 
energy gap. This type of behavior is demonstrated in Fig. [38^ ’ for one-dimensional lattices at small values 
of J/U. In higher dimensions, no exact results have been reported in the literature for the MI at T = 0. 
Nevertheless, high-order SCE jBn| as well as the quantum rotor approach supplemented by the Bogoliubov 
theory |3fi2| reveal exponential decay of Fa for small J/U in two and three dimensions. In Refs. [80113fi2] it 
was also shown that the correlations along the lattice diagonals are weaker than those parallel to the lattice 


axes in agreement with Eqs. (261l-(263) and numerical data in Fig. 38 d. 


In the SF phase, the two-point correlations decay according to a power law. In a one-dimensional 
lattice, the one-body density matrix does not show off-diagonal long-range order and its long-distance 
asymptotics [TZl I281| 


Fa{s) « 


(314) 


is determined by the Tomonaga-Luttinger parameter (2931. The prefactor A in the case of hard-core bosons 


coincides with the coefficient C in Eq. (2471 and for the ideal Bose gas A = (hi). This sort of behavior 
demonstrated in Fig. |3^a’) was studied in details in Refs. |318| using the TEBD algorithm. 

In higher dimensions. Fa takes constant values at large distances in the SF phase, which is a manifestation 
of the Bose-Einstein condensation. This was confirmed by QMC calculations for two-dimensional systems 
in discrete |36011363| and continuum |175| models. 

The finite-size effects can be well controlled by comparison with the SCE. This is demonstrated in 
Fig.HKb) for two-dimensional systems, where one can see that exact diagonalization for the lattice of 3 x 3 
sites overestimates the strength of correlations, while QMC data obtained for 10 x 10 sites show a very good 
agreement with the SCE. In three dimensions, QMC data for Fa{l) obtained for the lattice of 5^ sites are 
not affected by the finite-size effects. 

The quasi-momentum distribution P{k) in the ground state calculated for a one-dimensional lattice in 
the case of unit filling by exact diagonalization is shown in Fig. In general, it is periodic and even 
function of k which monotonically decreases in the interval k G [0,7r/a]. For small J/U, the distribution is 
rather fiat and this situation is well described by the perturbative expansion in J/U. With the increase of 
J/U, the maximum at fc = 0 grows and becomes very sharp for large J/U. Quantitatively, the form of the 


quasi-momentum distribution is characterized by visibility defined in analogy to Eq. (1091 as 


V = 


P(0) - P{7r/a) 
P{0) + P{TT/a) 
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(315) 






















Figure 37: (color online) Probabilities of the occupation numbers ni = 0 (red), 2 (blue), 3 (magenta) of bosons at a lattice 
site in the case of unit filling in one (a), two (b), and three (c) dimensions at zero temperature. Solid lines are the results of 
exact diagonalization for 14 sites in ID (a) and for 3x3 sites in 2D (b). Dots are QMC data for 10 x 10 sites in 2D (b) and 
5x5x5 sites in 3D (c), see also Fig. 1 in Ref. [2721 and Fig. 3 in Ref. I360l . Dashed lines in panels (a,b,c) are the results of 
the strong-coupling expansion up to the fourth order [Eq. ( |253| l]. (a’) Particle-number fluctuations in a one-dimensional lattice 
of 11 sites at different temperatures: k-QT/U = 0 (red), 0.16 (green), 0.2 (blue), 0.3 (magenta), 0.5 (brown). Red dashed line 
is the result of SCE at T = 0 up to {J/U)^^ |53| . 


It is a monotonic function of JjU which shows a linear dependence for small J/U and becomes close to 
unity near the quantum critical point. Quantum Monte Carlo calculations show that this feature remains 
preserved in higher dimensions and the measurement of visibility can be used in order to detect not only 
the quantum critical point but also the critical temperature of the transition from the superfluid into the 
normal phase |flh5] . 


9.7. Higher-order correlation functions 

In Figs.|40[|42[ we present exact numerical results for the particle-number and parity correlation functions 


Fn and defined by Eqs. (1311, (147). The correlations grow with the increase of J/U up to the 

maximal value and then decrease approaching the limit of the ideal Bose gas. Similarly to the one-body 
density matrix, the second-order correlation functions decrease exponentially with the distance in the MI 


phase which is demonstrated in Fig. 40 |a’) for in one dimension. This is consistent with the results of 
the SCE (see Eqs. (269l-(276l and Refs. |531150] ! and remains valid in higher dimensions. The second-order 
correlation functions decay also faster along the lattice diagonals than along the axes. 

The particle-number correlation function Fn is in general negative because due to the conservation law 
the increase of the particle number on one lattice site should be compensated by the corresponding decrease 
on the other site. In one dimension, the large-distance behavior of the particle-number correlation function 
in the SF phase is given by |281| 


Fn{s) « - 


1 


A{hif' cos {2'K{ht)s) 




((h^)s) 


2//Ct: 


(316) 


where the second term does not vanish only for incommensurate fillings {hf). From Eqs. (222), (241) it 
follows that the particle-number correlation function of hard-core bosons in the thermodynamic limit and 
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Figure 38: (color online) Correlation function Fa{s) at zero temperature in one (a,a’), two (b), and three (c) dimensions for 
unit filling, (a) Solid lines are obtained by the DMRG calculations in ID for the lattice of 128 sites and s = 1 (red), 2 (green), 
3 (blue), 4 (magenta). Dotted lines are the results of the strong-coupling expansion up to the 13th order m- (b) Solid lines 
are the results of exact diagonalization in 2D for the lattice of 3 x 3 sites and s = 1 (red), y/2 (green). Dots are QMC data for 
the lattice of 10 x 10 sites and s = 1. (c) Dots are the results of QMC calculations for the lattice of 5 x 5 x 5 sites and s = 1. 
Dashed lines in panels (a,b,c) show corresponding results of the strong-coupling expansion up to the third order in J/U, see 
Eqs. ( |261| l- |263| . (a’): Dependence of Fa on s in ID for J/U = 0.2,0.25,0.3,0.4,0.5 from bottom to the top obtained by the 
DMRG calculations, see also Refs. [2501131^ I364j . 



kir/a J/U 


Figure 39: (color online) (a) Quasi-momentum distribution at zero temperature in a one-dimensional chain of L = 14 sites 
with periodic boundary conditions in the case of unit filling for J/U = 0.05 (red), 0.1 (green). Dots - exact diagonalization, 
dashed lines - results of the strong-coupling expansion up to the third order in J/U, see Eqs. ( |265| l, ( |266[ |. (b) Visibility of the 
interference pattern for the same system obtained by exact diagonalization (solid line) and strong-coupling expansion (dashed 
line). 
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Figure 40: (color online) Particle-number correlation function Fn{s) at zero temperature in one (a,a’), two (b), and three (c) 
dimensions for unit filling, (a) Solid lines are obtained by exact diagonalization in ID for the lattice of 14 sites and s = 1 (red), 
2 (green), 3 (blue), 4 (magenta). Dotted lines are the results of the strong-coupling expansion up to the 14th order m- (b) 
Solid lines are the results of exact diagonalization in 2D for the lattice of 3 x 3 sites and s = 1 (red), y/2 (green). Dots 
are QMC data for the lattice of 10 x 10 sites and s = 1 (red), y/2 (green), (c) Dots are the results of QMC calculations 
for the lattice of 5 x 5 x 5 sites and s = 1 (red), y/2 (green). Dashed lines in panels (a,b,c) show corresponding results of 
the strong-coupling expansion up to the fourth order in J/U, see Eqs. ( |269[ |-( |27Tt . (a’): Dependence of \Fn\ on s in ID for 
J/U = 0.2,0.25,0.3,0.4,0.5 from bottom to the top obtained by DMRG calculations for 128 sites. 



kajiT 


Figure 41: (color online) Static structure factor at T = 0 in a one-dimensional homogeneous lattice of L = 14 sites with 
periodic boundary conditions calculated by exact diago nalization for unit filling and for J/U = 0.1. 0.2. 0.3. 0 .4. 0.5.1. 5.10 
(from bottom to the top). So{kq) is defined by Eq. (130i. Soikq) is obtained from So{kq) according to Eq. (129 i using Go{k) 
for Vq = 10 Lines are guides to the eye. 
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Figure 42: Parity correlation function at zero temperature in one (a), two (b), and three (c) dimensions for unit 

filling, (a) Solid lines are the results of exact diagonalization in ID for the lattice of 14 sites and s = 1 (red), 2 (green), 
3 (blue), 4 (magenta), (b) Solid lines are the results of exact diagonalization in 2D for the lattice of 3 x 3 sites and s = 1 (red), 
\/2 (green). Dots are QMC data for the lattice of 10 x 10 sites and s = 1 (red), y/2 (green), (c) Dots are the results of QMC 
calculations for the lattice of 5 x 5 x 5 sites and s = 1 (red), \/2 (green). Dashed lines in all panels show corresponding results 
of the strong-coupling expansion up to the fourth order in J/U, see Eqs. ( |274[ |-( [27^ . (For interpretation of the references to 
colour in this figure legend, the reader is referred to the web version of this article.) 


for ^ ^ I has exactly the form (3161 with /Ctl = 1 and A = l/(27r^). 

Static structure factor at T = 0 in a one-dimensional lattice in the case of unit filling is shown in Fig. 

If J/U is sufficiently small, S'o(fc) as well as 5o(fc) are quadratic functions of k for small k. This can be 
clearly seen in spite of rather coarse discretization in the momentum space . For larger J/U, So{k) and 
So{k) become linear functions of k for small k in agreement with Eq. (1451. This indicates an advent of 


the sound mode making a dominant contribution into the structure factor. In this regime, with the aid of 


Eqs. (238), (293) we can rewrite Eq. (145) in the form 


lim So{k) = 


k 


k^O 


2fcF^TL 


(317) 


where kp is defined by Eq. (235), which shows that the behavior of the structure factor in the limit of small 


k gives an access to the Tomonaga-Luttinger parameter. Using the condition /Ctl = 1/2 (see section 8.1), 


one can determine the critical value oi J/U for the transition from the superfiuid into the Mott insulator. 
Calculations of So{k) by exact diagonalization for chains up to L = 14 sites give {J/U)c « 0.28 |366| in 
perfect agreement with other calculations based on the exact diagonalization listed in Table 

The parity correlations are of the same order of magnitude as the particle-number correlations. However, 
the former possess a more narrow maximum near the critical point. In one dimension, this maximum is 
below {J/U)c and captured by the SCE of the 4th order for the nearest-neighboring sites. In two and three 
dimensions, the maximum is above {J/U)c and, therefore, out of reach of the perturbation theory in J/U. 
Exact results for one and two dimensions presented in Fig. [4^a,b) are very similar to that of Ref. |194| . 
where one can also find a comparison with the experimental data obtained at finite temperatures in the 
presence of the harmonic confinement. Positive values of Fg are due to the fact that a change of the particle 
number on site R compensated by an opposite change on site R leads to the same change of the parity on 
both sites. 


76 




















10. Mean-field theory 

Mean-field theory of lattice bosons relies on the concept of spontaneous breaking of U(l) symmetry. 
According to the Mermin-Wagner-Hohenberg theorem |36711368] it is valid in two-dimensional systems at 
zero temperature and in higher dimensions at arbitrary temperature. It gives an exact solution in the limit 
d ^ oo |116l 111711369| . implying that J —>■ 0 such that dJ is finite, as well as in the case J = 0 and arbitrary 
d. 


10.1. Decoupling approximation 

The mean-field theory is based on the assumption that the second-order fluctuations of the bosonic 
creation and annihilation operators are negligible, i.e., — (“la)) ~ ^or li ^ I 2 . This leads 

to the decoupling approximation of the hopping terms |471 IT51 1115L1339] 


A 




t 


Then the Bose-Hubbard Hamiltonian (541 becomes a sum of local operators 

d 

Hbh ~ ~ [“1 + A*+eJ - V'l + h-C. 

1 I/=l 


u 


^afaj'aia, 


(318) 


(319) 


where i/'i = (“i) are c-numbers. The neighboring lattice sites in the Hamiltonian (3191 are coupled only 
via the expectation values of the creation and annihilation operators. Thus, the mean-field theory neglects 
quantum correlations between different sites and the states of the Hamiltonian (319) are tensor products of 
the local states 


l®) = (^|si). |si) = 


(320) 


n—0 


which is equivalent to the Gutzwiller ansatz |75l 136911370] . Here |n)i is the Fock state with n atoms at site 1. 
Normalization of the |si) imposes 

00 

n—0 


As it follows from the form of the state (320), the Gutzwiller approximation neglects quantum correlations 
between different lattice sites but takes into account on-site quantum fluctuations, provided that |si) is not 
a single Fock state. 

The mean number of condensed atoms on a lattice site 1 in this approximation is given by I'l/'ipj where 


V'l = (ai) = 


iQn V ^ 


(321) 


is the condensate order parameter. One can easily show that IV'iP cannot be larger than the mean occupation 
number 


(^i> = 


n \Cln\ 


(322) 


Minimization of the functional 


- cindtCiJ - (H^i) + ^iN 

n—0 
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leads to the system of Gutzwiller equations (GE) |37111372| : 


ih 


dc 


In 


dt 


= 


Qn' 


^nn 


n'—O 


u 


—n(n — 1) — 


( 323 ) 


Sn',n - JVTT/Sn',n+1 ^ (V'f+e^ + V^r-e^) “ JV^Sn,n' + l ^ (V'l+e, + '01-e„) ■ 


Note that these equations are invariant under transformation ci„ —>■ (—l)"ci„. 

Although the mean-field Hamiltonian (319) does not satisfy all fundamental commutation relations of 
the original Bose-Hubbard Hamiltonian (54), the Gutzwiller approximation can be considered as conserving 
because the expectation values of the total number of particles, total energy, and the quasi-momentum 
remain constant in time. Equations (323) allow to study the properties of the ground state as well as the 
dynamics of excitations. 


In the mean-field approximation, the one-body density matrix takes the form 

(aqab) = V'qV'b +<^1112 “ l^ul^) • 


(324) 


In a homogeneous lattice, and the largest eigenvalue of the one-body density matrix in the ther¬ 
modynamic limit is Nq = L‘^ Therefore, the condensate fraction /c = /(ni) which coincides 

with the superfluid stiffness /g. This implies that Eq. (312) is valid in this case. In an inhomogeneous lattice, 
one cannot write in general an explicit expression for Nq and it has to be determined numerically. This can 
be efficiently done by the iteration procedure |373| 


N, 


(®+i) _ Arid 




IV'i 


- {n\) + 1^11 


(325) 


obtained from a rearrangement of the eigenvalue equation ( |99[ ). 

As it was shown in Ref. IM], in the weakly interacting regime, C//J <C 1, Eq. ( |323| ) can be transformed 
into the discrete Gross-Pitaevskii equation (DGPE) 




dt 


(326) 


i/=i 


assuming that |si) in Eq. (320) are Glauber coherent states (see also Ref. |375| 1. Eq. (326) describes a pure 
Bose-Einstein condensate in a discrete lattice model, which implies (ni) ~ |r/'i| . 

10.2. Ground state 


In the homogeneous lattice, the ground state of the Hamiltonian (319) is described by the coefficients 
ci„ that do not depend on the site index 1. This corresponds to the stationary solution of Eqs. (323) of the 
form 

cin{t) = exp {-iujot) , (327) 


tvjjQ = —idJ 


v.'»> +5: 


n—0 


—n[n — 1) — fin 




(328) 


The explicit form of Cn'’ depends on the Ailing factor (ni) and the ratio J/U. If the latter is less than the 
critical value determined as HHl 


2d(J/C/)e 


(no - fi/U){fi/U - no -I- 1) 
1 + fi/U 


(329) 
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Figure 43: (color online) Probabilities of the occupation numbers n = 0 (red), 1 (green), 2 (blue), 3 (magenta) in the ground 
state for the filling factor {hi) = 1. Solid lines are exact mean-field results and dashed lines show the asymptotic values 
in the limit U ^ 0 according to Eq. ( |334| l. 


for no — 1 < fi/U <no, the solution is given by 


^(0) _ X 

■"r, - Oi 


n,no 5 


(330) 


where ng is the smallest integer greater than ^/U. In this case the superfluid order parameter = ipx 


deflned by Eq. (3211 vanishes and we have the MI phase with exactly ng particles at each lattice site. 


Eq. (329 \ determines the dependence of the critical ratio of J/U on the chemical potential. It can be also 


inverted to determine the dependence of the critical chemical potential on J/U, which gives two solutions 


u±(no) 1 /. 2dJ\ , 1 /, 

%^ = n„--h+—\±-Jl-(4.n„ + 2)u 


2dJ ^ f2dJ 


V u 


(331) 


that are real, provided that 


2d{j/u) < 2d(j/t/)r" = 


For J/U = (J/U)^^^, the two solutions merge into one 


^^±{no) 

U 


= (^)^ = Vn-o{no + 1 ) - 1 


(332) 


(333) 


which describes the tips of the Ml-lobes on the phase diagram. 


If we expand Eq. (3311 up to the third order in 2dJ/U, we immediately observe that it coincides with the 
results of the strong-coupling expansion given by Eqs. ( 281| ), (2821 in the limit Z = 2d ^ oo. This confirms 
that the mean-field theory becomes e xact in infinite dimensions. 


If J/U exceeds the critical value (3291, Cn'^ has a broad distribution leading to the fact that the order 


parameter does not vanish which corresponds to the SF phase. In this regime, exact analytical solution 
for Cn'^ can be obtained only in the limit (7 = 0 and has the form 


= exp I - 


,(..177^7;. he. 

/ vn! 


= (hi) , /r = —2dJ . 


(334) 


Eq. (334) describes the coherent state and the corresponding particle-number distribution |c: 
with the exact result (160) for the ideal gas. 


,coh I 


coincides 


For finite J and U above the critical value (|329|) exact coefficients Cn^ can be obtained by numerical 


diagonalization and the results for (hi) = 1 are shown in Figs. As small values of J/U, the 
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Figure 44; Mean-field phase diagram, (a) Green area is the MI region with uq = 1 described by Eqs. ( |329[ l, ( |331[ l. Thick 
green line shows the dependence fi{J) in the superfluid phase for (n) = 1. Other lines show also /i(J) for non-integer fillings 
which are multiples of 0.1 (see also Fig. 13 in Ref. |72| V In the gray areas, where > 0 at constant (h) the superfluidity 

has a hole character. In the rest part of the diagram, we have a particle SF. (b) Green areas show the first three MI zones 
(no = 1,2,3). The lines of constant (hi) are labeled by the corresponding filling factors. Gray areas are the regions of the 
hole superfluidity. In the rest part of the diagram, we have a particle SF. (For interpretation of the references to colour in this 
figure legend, the reader is referred to the web version of this article.) 


probabilities of having n = 0,2 particles nearly coincide and the probabilities of the occupation numbers 


larger than 3 are negligible. The comparison with the coherent-state distribution (334) shows that rather 


large values of J/U are needed in order to approach this asymptotic limit with a good accuracy. 

The quasi-momentum distribution (73) of the atoms in the ground state is given by (see, e.g., Ref. |376] ) 


PM = 


W(k) 


(^i) 


f- 

V27r 




V' 


( 0 ) 




(335) 


where the lattice is assumed to be infinite. The whole quasi-momentum distribution is confined within the 


momentum distribution 


W{k) 


of a single site. The first term in the brackets, {n\) — , describes 

incoherent part of the system, and the J-peaks at the vectors of the reciprocal lattice k = are clear 
signature of the Bose-Einstein condensation. Similar structure of the quasi-momentum distribution, although 
with somewhat smeared maxima instead of sharp peaks, was observed in the experiments 0131211 HI El 

HTHllIBniEZZl. 

As discussed in Refs. |3781 [379], near the phase boundary one has to distinguish between particle and 
hole superfluidity. For the hole SF, the function ^{J) at constant filling factor (n) has a positive derivative 
This is only possible for fillings uq — 0.5 < (n) < no as is demonstrated in Fig. 44 showing the 


corresponding hole SF regions. For the particle SF, on the other hand, < 0. Consequently, far away 

from the phase boundary, superfluidity has always a particle character. With the increase of the filling 
factor, the size of the regions of the hole superfluidity decreases together with the size of the Ml-lobes. As 
we will see in section [10.7| the difference between the particle and hole superfluidity plays an essential role 
for the character of the topological modes. 


10.3. Excitations 

We consider small perturbation of the ground state 


C\n{t) = + C^in^t) -k . . . exp (-IWqO , 


(336) 
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where 


CjW (t) = ■ (337) 

Substituting this expression into GE and keeping only linear terms with respect to Ukn and Ukn, we obtain 
the system of linear equations EHO], 


feuk 


Mk 

Wk 


73k 

—73k —^k 


Uk 

hk 


(338) 


where ifk and Uk are infinite-dimensional vectors with the components Ukn and Ukn ("n = 0,1,...), respec¬ 
tively. Matrix elements of and have the form 


— — Sn\n+1 + V^^n,n' + l^ + 


u 


n(n — 1) — /rn — fvjjQ 


^n' ,n 


-rynn 


— JVi 
= -Jk 


Vn + l Vn' + 1 1 -f 


n+l *-71' + ! 


.(0) JO) 


^n—1 '^n' —1 


Vn+l Vn' -f Vn Vn' + 1 c® ^ 


/ + ! 


where Jk = ‘^dJ — Ck with Ck being the energy of a free particle (1401. This system is valid for both phases 


and generalizes the Bogoliubov-de Gennes (BdG) equations previously derived for coherent states |3811I382] . 
The dependence on the vector k is determined by the variable 


i. 

I 1 • 2 \ 


(339) 


which varies from 0 to 1. For small |k|, x « \'k\a/{2'/d). 

The energy increase due to the perturbation is given by m 

AE = fiwk (|uk|^ - |uk|^) . 


(340) 


Formally, Eqs. (3381 have solutions with positive and negative energies i/iWk, which are equivalent because 
Eqs. (3771, (3401 are invariant under the transformation Wk —t —Wk, k —>■ —k, ftk —t v^, t Uk, so that 


only solutions with positive energies will be considered in the following. The eigenvectors are chosen to 
follow the orthonormality relations 


i.A' • Uk,A - 'f^,A' • Wk.A = I^A.A' ■ 


Perturbation (3771 creates plane waves of the order parameter '0i(t) = where 

= ^J^ghk x.-a;kt) y*g-i(k xi-u;kt) ^ 
oo 

TJk = (Mk,n+1 Vn -I- 1 -I- Uk,n-1 Vn) , 


n—Q 

oo 


Vk = (Mk,n-lA/» + fk,n-nV?^+ l) ■ 

n—0 

The perturbations for the total density and the condensate density are given by 

(n,)(t) = (hi)(°)+ [Mke*(*^'"'-‘"‘'*)+c.c.] , 

OO 

Ale — ^ ^ ^ (Ukn T "^kn) ^ 7 


(341) 


( 342 ) 


n—O 
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and 


\Mt)\ = 




(0) 




(343) 


Sk = {Hk + Vk) . 

In what follows we consider the properties of the excitations in the MI and SF phases. 


10.3.1. Mott insulator 


For the MI phase, the coefficients have a simple analytical form (330). The eigenvalue problem 


for the infinite-dimensional matrices (3381 reduces to the diagonalization of two 2 x 2-matrices that couple 


Uk,no-i to Wk,no-i-i a^d Uk,Tto+i to Wk,no-i> respectively. The lowest-energy excitation spectrum consists of 
two branches 


f^k± = 


1 


^C/2-JkC/(4no + 2) +J2± 


Ulno- 


Jk 

Y 


0 


(344) 


The same result was obtained using the Hubbard-Stratonovich transformation ma and within the Schwinger- 
boson approach |383| . 


These two branches are shown in Fig. 45 and display a gap. The sign in front of /i in Eq. (344) is different 


for the two modes. Therefore, the solutions labeled by ’-)-’ and ’ correspond to the situations, when one 


particle is added into the system and removed, respectively. Hence, the eigenmodes described by Eq. (344) 
are called particle and hole excitations |251| . If the total number of particles is conserved, the two kinds of 
excitations can be created only in pairs and the corresponding energies are added. 

Nonvanishing coefficients for the two modes are given by 


,(+) 

k,no —1 


^k,no + l 


1 2 


- 1 = 


k,no + l 


n 2 


'^k,no —1 


n 2 


-1 


(345) 


U - 2 Jk (no + i) 


yC/2-4JkC/(no + i) + J, 


- 1 


According to Eqs. (343), (342), no density wave is created in the two modes. However, the order parame¬ 


ter (341) does not vanish and takes the form 




‘'1+ 

',(1) 


YYit) = VYe 




(346) 


In the complex plane (Re('!/'),Im(i/))), this corresponds to the motion on the circles with radii |f^k+| and 
around tp = 0. 

Other solutions of Eq. (338[) are independent of k with the energies 


hujx = ^ [A(A - 1) - no(no - 1)] - - no) 


(347) 


They are denoted by A which are non-negative integers different from no, no ± 1. Since no is greater than 
/r/C/, the excitation energies are always positive. The eigenvectors of these modes have the form UknA = ^n,A> 
"^^knA = 0, and the amplitudes of all the waves defined by Eqs. (341), (342), (343) vanish. 

In Fig. [4^b) we compare the energies of the particle and hole excitations 


ek± = ± (/iWki ± m) 


(348) 


determined by Eq. (344) with the analogous quantities 


^kO 


- e: 


00 


^k± — ± {^noL‘^±l ^noL'^J 
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Figure 45: (color online) Excitation spectrum of the Mott insulator, (a) First three branches (1)^ ^•^k+ (2) and hux=Q (3) 

of the excitations spectrum of the MI for ^/U = 1.2 and 2dJ/U = 0.05, which corresponds to uq = 2 [Eqs. piil i, PStI i]. 
(Reproduced from Ref. ES], (c)2011 American Physical Society.) (b) Energies of particle and hole excitations p48| l calculated 
for J/U = 0.08 in one-dimensional lattice in the case of unit filling. Solid lines — mean-field theory [Eq. p44[ |], dashed lines - 
strong-coupling expansion [Eqs. (|283[l, (|284[|], dots - exact diagonalization for N = L = 13. 


calculated by exact diagonalization and strong-coupling expansion for uq = 1 and d = 1. We see that the 
strong-coupling expansion of the 6th order is in perfect agreement with the exact numerical results. The 
mean-field theory based on the Gutzwiller ansatz gives qualitatively correct predictions but underestimates 
the gap for the particle-hole excitations. 

The boundary between the SF and MI phases is determined from the disappearance of the gap in the 
excitation spectrum, i.e., when wq- = 0. Under this condition, we recover the critical ratio (3291. For 
J/U > {J/U)c, the lowest frequency wq- in Eq. (3441 becomes negative leading to a negative expression for 


Eq. (3401, so that the Mott-phase solution (3301 does not correspond to the ground state anymore. 


The excitation spectrum has peculiar features on the boundary between the MI and SF. For { J/U)c = 
the excitation energies (3441 can be rewritten as 


/iCJki = 


\/no(no + 1) Uck + ^ 


-I 1/2 


-I 


(350) 


For small |k|, the two branches are degenerate and have linear dependence Wk± = Cs’^|k| with the sound 
velocity 


c“P = 


U -v/ncT+T- yTid 


[noino + 1 )] 


1/4 


(351) 


hy/d 

expressed in units of the number of sites per second. For other points on the boundary, i.e., {J/U)c < 
{ J/U)y’^^, no degeneracy appears and the sound velocity vanishes leading to the quadratic dispersion a;k± 
for small Ikl. 


From the expression for the excitation energies (3441 one can obtain the values of the mean-field critical 


exponents for the MI-SF phase transition. According to the scaling theory m, the energy gap should be 
proportional to \t — where t is a control parameter which approaches its critical value tc, ^ is the 

dynamical critical exponent and v is the critical exponent of the correlation length. 

First we consider the transition at fixed J/U under variation of the particle number. In this case, the 
role of the control parameter is played by fi. Since the energies fiwki in Eq. ( 344[ ) are linear functions of /i, 
we get zv = 1 which is consistent with the fact that z = 2 and v = \/2. 


Now we consider the transition at fixed particle number controlled by the ratio J/U . Using Eqs. (3441, (332), 


we can write the expression for the total energy required to create one particle-hole excitation at k = 0 in 
the form 

1/2 / _ 1/2 


hujQ— -fi hujo-\- — u 


4 ^/ no (no -fi 1 ) -fi Jc — J (jc ~ 


(352) 


where J = 2dJ/U. When J is close to Jc, the excitation energy is well approximated by the lowest-order 
term 

hujo- -fi hujo+ « 2U [no(no -fi 1)]^'^"^ ( Jc — j) , (353) 
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Figure 46: First three branches \ (A = 1, 2, 3) of the excitations spectrum of the SF for fi/U = 1.2 and 2dJJV = 0.15. The 
straight dashed line represents the linear approximation with the sound velocity | |356| |. (Reproduced from Ref. |374| . ©2011 
American Physical Society.) 



Figure 47: Flatness parameter ||355^ for ^/U = 1.2 and 2dJ/U = 0.15 (a), 1 (b). 


which gives zi/ = 1/2. This is consistent with the result z = 1, = 1/2 for the classical XY-model in the 

dimensions larger than three. 


10.3.2. Superfluid 

In the SF phase, the eigenvalue problem (3381 is solved using the numerical values of c\ 
and p./U. The energies of the lowest-energy excitations are shown in Fig. 46 
consists of several branches which form a band structure. 


for each J/U 
The excitation spectrum 


We note that only the first two lowest-energy 
branches have a strong dependence on k. In the complex plane (Re('!/'), Im(^)) different modes correspond 
to the motion around 0 on the ellipses with the axes 

&kA = l^kA + I , 6L = l^^kA - ViJy I . (354) 

In order to distinguish between different modes, it is convenient to introduce a flatness parameter |384| 

/kA = (&kA ~ ^kA)/(^kA + ^kv) : l/kA| < 1 ■ (355) 

In contrast to the MI, the lowest branch has no gap. It is a Goldstone mode that appears due to the 


spontaneous breaking of the U(l) symmetry. The flatness parameter (3551 for this mode is negative which 
is interpreted as phase-like oscillations |38411385] . As is shown in Fig. 48 the amplitude of the total-density 
wave is larger than the amplitude of the condensate-density wave for this mode. A value greater than unity 
for the ratio Aki/Bki means that the condensed part and the normal part oscillate in phase. 

The lowest-energy branch has a linear form Wk,i = |k| for small k with the sound velocity given by |374| 


c° = 


iP 


(0) 


/h, 


(356) 


where k is the compressibility (82 b This result proves that the Gutzwiller approximation is gapless. 
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Figure 48; for = 1-2 and 2dJjU = 0.15 (a), 1 

Society.) 



(b). (Reproduced from Ref. | 374 |. ©2011 American Physical 



Figure 49: (a) Sound velocity calculated numerically from Eq. ( |356| |. Note the discontinuities at the points {fi/U)c] 

described by Eq. ( |351[ l. (Reproduced from Ref. [3741 . ©2011 American Physical Society.) (b) Comparison of the sound velocity 
calculated numerically from Eq. l|356|l (solid lines) with the analytical expression ||358|l (dashed lines) for (hi) = 0.5,1. 


Fig. 49 shows the dependence of the sonnd velocity on and J calcnlated nnmerically nsing Eq. (3561. If 


we approach the boundary of the MI, the sound velocity goes to zero everywhere except the tips of the lobes, 


where it is perfectly described by Eq. (3511. This behavior can be understood considering the properties of 
and K. If we approach the SF-MI transition from the SF part of the phase diagram, the order parameter 
tends always continuously to zero. The compressibility k reaches a finite value at every point of the 
boundary except the tips of the MI-lobes where it tends continuously to zero such that the ratio 
is finite. Therefore, the sound velocity vanishes at any point of the phase boundary except the tips of the 
lobes 071 [ 73 . 

For a weakly interacting gas {U <C J), « (hi) and n « 1/(7. In this limit, we recover the 

Bogoliubov dispersion relation 

huj^ = y/ck (ek + 2(7(hi)) (357) 

and the expression for the sound velocity |38611387] , 


cf = y /2 J(7(hi )/h . 


(358) 


A comparison with the exact numerical values calculated according to Eq. (3561 shows that the approxima¬ 


tion (3581 overestimates the sound velocity and predicts completely different behavior at small tunneling 
rates and integer fillings [see Fig. @b)]. 

In the opposite limit (J ^ U), the sound velocity is given by |374| 


c° = 2J(no -f l)\/2d((hi) - n.o)(no + 1- (hi))/h . 


(359) 


It vanishes at (hi) = no,no + 1 and takes maximal values at (hi) = no -I- 1/2. This qualitative behavior is 
the same as in the case of hard-core bosons in ID, where the sound velocity is given by Eq. (233). 
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Experimentally the speed of sound can be measured with the aid of an external potential which creates 
a density perturbation of the gas |8881 1889] . Corresponding numerical simulations of the sound waves 
propagation for the lattice Bose gas were performed within the framework of the mean-field theory on 
the basis of the DGPE |386| as well as GE |374| and show perfect agreement with Eqs. (358) and (356), 


respectively. Exact numerical simulations for soft-core bosons in ID were also done making use of the DMRG 
method |355| . 

Another possibility is to employ the Bragg spectroscopy which gives an access to the momentum-resolved 
excitation spectrum. The experimental data for the lowest excitation branch of the superfluid obtained in 
Ref. |40| agree qualitatively very well with Eq. (357). However, the energies for small momenta appear 


to be slightly underestimated and for high momenta overestimated which might be an indication that the 
calculations beyond the Bogoliubov theory are needed in this regime. 

Higher modes (A > 2) have gaps that grow with the increase of J. For the second mode (A = 2), the 


flatness parameter (355) is positive (see Fig. 47). This type of excitations corresponds to the amplitude-like 


oscillations of the order parameter which is called ‘Higgs’ amplitude mode |3831138511390] . The amplitude 
of the total-density wave is much less than that of the condensate-density wave (see Fig. 48) meaning that 
the oscillations of the condensate and normal components are out-of-phase. The properties of this mode 
have been studied in details in the context of the Bose-Hubbard model within the framework of the mean- 
field theory |7S 1 1771 IHTil IMnHMK] and in quantum Monte Garlo calculations |39()H398] . It was 

experimentally observed with the aid of Bragg spectroscopy in the non-linear regime |390| and with lattice 
modulation in the linear-response regime |5M] . However, the DMRG calculations of the dynamic structure 
factor in one dimension did not reveal any distinct gapped modes in the superfluid phase |2501139^ . 


10.4- Bragg scattering 

Within the Gutzwiller approximation, the susceptibility in the lattice version of Eq. (116) can be written 
in the form |374| 


X(k,w) = 


XkAWkA 


{oj + iO) — wj 


(360) 


kA 


where A denotes various excitation branches discussed in the previous section associated to the eigenvalues 
WkA and the corresponding amplitude of the Bragg scattering is determined by the amplitude of the density 


wave (342) as 


XkA = lAkAT 


(361) 


The dependences of Xk,A on the variable x defined by Eq. (339) for the excitation branches with A = 1,2,3 
in the SF phase are shown in Fig. For the chosen values of parameters, only the two lowest branches 
display noticeable amplitudes. In the long-wavelength limit, only the lowest mode has a nonvanishing 
amplitude 


Xk,i 


= ic°|k| . 


(362) 


Similar results have been also obtained in Ref. |383| . However, the calculations in Ref. |383| are valid only 
close to the boundaries of the MI-SF phase transition because the occupation numbers n in Eq. (320) were 


restricted to n = no, uq ± 1. In Fig. [^instead, we see that the amplitude for the third excitation branch as 
well as for the second one can become significant at certain densities. 


As in the case of a Bose gas in continuum, the sum rules (136), (139), (141) allow us to deduce the static 
structure factor. We And indeed 


Soik)=J2 


Xk,A 


= 2Cs|k| 


(363) 


This result shows an interesting feature of the sum-rule approach. Starting from the lowest-order Gutzwiller 
approximation that does not contain any correlation, the two-point correlation function is determined as a 
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Figure 50: Transition amplitudes Xk.A associated with the transition frequency ojk for the lowest excit ation branches 

(A = 1,2,3) and for ^ijU = 1.2 and 2dJ/U = 0.15. The dashed line corresponds to the approximation | |362^ . (Reproduced 
from Ref. 13741 . ©2011 American Physical Society.) 
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Figure 51: (color online) Transition amplitudes Xk,A versus the density for the first excitation branches (A = 1,2,3) and for 
X = 1 and 2dJ jlj = 0.2 (a), 0.3 (b). Dashed lines show the static structure factor 5o(k). (Reproduced from Ref. |374| . ©2011 
American Physical Society.) 




next-order contribution. Similarly, starting from the time-dependent DGPE, an analogous procedure has 
been successfully used to recover the static structure factor predicted from the Bogoliubov theory |1711I400] . 

In the MI phase, the Gutzwiller approximation does not allow us to observe any branches since Xk,A = 0. 
No Bragg response is possible, although the excitations exist in the mean-field approach. In order to allow 
a nonvanishing response, correlations between different sites should be included, which goes beyond the 
Gutzwiller approximation |383L I401| . In such a description, excitations in the Bragg process are created 
as particle-hole pairs I3H3 uni IMI- However, the latter appears to be of the second order in the inverse 
of coordination number Z = 2d |40111403H405| and, therefore, is not taken into account by the standard 
Gutzwiller ansatz. 


10.5. One-particle Green’s function 

In the context of the Gutzwiller approximation, the one-particle Green’s function can be determined as 
a linear response to the perturbation |374| 

H’{t) = +h.c. , (364) 

1 


which explicitly breaks the U(l) symmetry. This induces fluctuations of the order parameter. 


/ '01 — 

V 01* - (0^°^)* 


= G(k,cc). 


-iCk-xi-ujt) 


(365) 


The proportionality term is the one-particle 2x2 matrix Green’s function with the general expression 


0k,a;)=^ 

A 


^k,A 


uj -\-i0 — Wk.A 


(366) 
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where the matrix transition amplitude is defined as 


La = 






Wk,AVk,A 

|Vk,A|' 


(367) 


and the functions 7/k,A) Vk.A are determined by Eq. (341). 

In the MI phase, matrix g is diagonal and according to Eqs. (|341|), (|345|) its nonvanishing entries are 

given by 


=k,A 


|Vk-|^ = ^ 


(2no + 1)U - Jk 


yC/2-4JkC/(no + i) + J, 


- 1 


|Z7k+|" = |Vk_|" + l 


(368) 


Thus, the time-dependent Gutzwiller approach allows to recover the results previously established in the 
context of quantum field theory The spectral weight |Vk-|^ of the hole branc h yields the quasi¬ 

momentum distribution mm- For J = 0, |Vk-| = no in agreement with Eq. ( |335[ ), but for any finite 
J < Jc it has maxima at the vectors of the reciprocal lattice k = which grow as we approach the critical 
point. 


10.6. Finite-temperature phase diagram 

At finite temperature, the expectation values of the operators are calculated for the thermal state and 
the definition of the order parameter ipi should be generalized as 


V-i = 2^1 "(a*) 


i^s 


(369) 


where '0s,’s are defined by Eq. (321) for a particular state |si) with the eigenvalue of the local mean-field 


Hamiltonian (319) in the grand-canonical ensemble and Zi{p) is the corresponding partition function. In 


the homogeneous case, all the 0i are equal to each other. The region in the (/r, J) plane with vanishing 0 
is described by the equation |406l 14071 


jj ^0 \F) 


(l + s/Cf)exp(-to!) 

{n- fi/U){g/U -n-^1) 


OO 


< 1 


(370) 


where the partition function Zq{p) and the energies En are exactly the same as in the absence of hopping 


[see Eq. (179)]. In the limit of vanishing temperature, Eq. (370) reduces to (329). 


The boundaries of the superfiuid-insulator transition described by Eq. (370) are shown in Fig. 52 for 


different temperatures. The size of the insulating regions grows with temperature and the topology is the 
same as in QMC simulations [see Fig. [2^a)[. The insulating region is divided into two parts, MI and normal 
gas, separated by the crossover lines determined from the condition that the compressibility is fixed by a 
small arbitrary number. In the mean-field approximation, the latter is given by 


= ((^f) - /(kBT) 


(371) 


which follows from Eq. (83) if we neglect correlations of the particle numbers at different sites. Since in the 


mean-field theory the properties of the insulator do not depend on J, the crossover lines are parallel to the 
J/U axis in contrast to the QMC calculations. 

Eq. ( 370| ) allows also to determine the critical temperature of the superfiuid-insulator transition. For 
a given filling (hi), the chemical potential can be eliminated with the aid of Eqs. (180), (181) and the result 


for Tc is shown in Fig. 53 In the case of unit filling, the dependence of on U is qualitatively similar to that 


obtained in QMC calculations in three dimensions (see Fig. 29). However, the mean-field theory predicts 
larger values of compared to QMC and cannot reproduce correctly the ideal-gas limit (k^T^/J = 5.591). 
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Figure 52: (color online) Finite-temperature mean-field phase diagram. The boundaries of the superfluid-insulator tran¬ 
sition | |370| l are shown by solid lines for k^T/U = 0 (blue), 0.05 (green), 0.1 (red). Green stripes are the MI regions for 
k^T/U = 0.05 obtained from the requirement nil < 0.01. Outside of the stripes within the insulating phase we have a normal 
gas. 
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Figure 53: Critical temperatu re of t he su pe rfluid -insulator transition for the filling factor (fii) = 1 (a), 0.75 ( b). S olid lines 
are numerical solutions of Eqs. | |370|l, l |180| l, Dashed line in panel (a) is approximate analytical solution l |372| l obtained 

within the slave-boson approach 13401 . Dotted line in panel (a) is the crossover line between the MI and normal gas determined 
from the condition A = hLJo+ + fiuio— = k^T^, where huJo± are the energies of particle and hole excitations ^344[ |, see also 
Refs. [ 339 II 34 II . 


Other versions of the mean-field theory yield somewhat different results for the critical temperature. For 
instance, if the occupation numbers are first restricted to ng, ug ± 1 and then the mean-field approximation 
is applied, the critical temperature becomes lower, although the behavior of Tc(U) remains qualitatively 
the same m- Different form of Tc(U) was obtained in the slave-boson approach, where Tc{U) has a 
maximum fMT] . If the occupation numbers are again restricted to ng, ng ± 1, the latter approach gives an 
analytical expression |340| 



U-8 (2ng -b 1) 


U -h 2no + 1 


2dJ 2 

C-2(2ng-bl) 


U + A (2no + 1) 



(372) 


This result is also shown in Fig. 53 a) for comparison. 


10.7. Quantum solitons 


Superfiuid phase of the ultracold bosons in optical lattices far from the transition into 
described by the DGPE (3261 if the lattice is deep or by its continuum counterpart m 


the MI is well 
in the case of 
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shallow lattices. One of the remarkable features of these equations is that they allow soliton solutions in 
analogy to nonlinear optics |408| . This has triggered theoretical interest in discrete (lattice) solitons in the 
context of ultracold atoms |40911410] . and has led to the seminal observations of gap solitons, i.e., lattice 
solitons with repulsive interactions, but with an appropriate dispersion management |411| . 

While most of the studies of solitons were concentrated on their classical aspects, more recently, consid¬ 
erable interest has been devoted to the effect of thermal noise |410L 1412] , quantum properties of solitons, 
and the role of quantum fluctuations |41fl| . The latter may cause Ailing up of the dark soliton core in 
the quantum detection process, as was shown using the Bogoliubov-de Gennes equations |414| . The same 
method was also employed to study the stability of solitons |3811141511416] . excitations caused by the trap 
opening m, and entanglement generation in collisions of two bright solitons [418] . A noisy version of the 
standing bright solitons was studied using the exact diagonalization and quantum Monte Carlo method |419| . 
Bright solitons in ID were considered in Ref. mzi, where exact Lieb-Liniger solutions were used to calculate 
the internal correlation function of the particles positions. Making use of the DGPE, and the time-evolving 
block decimation algorithm |420| it was demonstrated that quantum effects cause the soliton to All in, and 
that soliton collisions become inelastic |421H423] . In the next section, we consider the properties of the 
discrete dark solitons near the SF-MI transition within the framework of the Gutzwiller mean-field theory. 


10.7.1. Standing modes 

In the present section, we consider low-energy excited states, where the coefficients ci„ as well as the 
order parameters ipi depend only on one spatial direction v. Without loss of generality, we can assume that 
this is i/ = 1. Then 

, if = 1, 

V'; , if > 1. 


V'l±e, = 


(373) 


We are interested in the stationary solutions of Eqs. (3231 


Cln{t) 

huji 


exp (-iw/t) , 


= -2J 


+ 2(d - l)^'^^ + 


/,( 0 ) 


(0) 


/,( 0 ) 


n—0 


u 


,(n — 1) — /in 


,(o) 


(374) 


where is determined by the coefficients according to Eq. (3211. We require that the condensate 
order parameter is an antisymmetric function with respect to the middle point of the lattice /q. These 
are the kink states which can be treated as standing dark solitons. In contrast to the ground state discussed 
in Sec. 110.2] all the quantities which describe the solitons are labeled by the site index 1. 

In general, one has to distinguish between the two cases: when the middle point Iq is on the lattice site 
(on-site modes) and in the middle of two neighboring sites (off-site modes). The two modes have different 
energies, and the difference defines the Peierls-Nabarro barrier |42411425] , which may affect the mobility of 
solitons. In addition, the stability of the on-site and off-site modes can in general be different, as it will be 
shown in the next section. 

We consider first the SF phase. Typical behavior of the kink modes with the lowest energy is displayed 


in Fig. 54 |379| . The mean occupation numbers {fii) calculated according to Eq. (3221 are shown in (a) 


and (c), while (b) and (d) give the associated defined by Eq. (|321|. The individual curves correspond 


to different tunneling rates J. Far from the middle point of the lattice, (n;) as well as tend to the 
same values as in the ground state. Near the middle point, on the other hand, they have nontrivial position 
dependence. 

For the considered chemical potential, /r/D = 1.2, the MI-SF transition occurs according to Eq. ( 329| ) 
at 2d{J/U)c ~ 0.0727. Much above this value, (n/) has only one extremum which is a global minimum. 
It is doubly degenerate in the case of the off-site modes [Fig. [54)(a), curves (i)-(iii); Fig. (He), curve (i)]. 
Expectedly, these solutions reproduce the well-known standing soliton of the DGPE |424L 1425] . For smaller 
values of J, when we come closer to the phase boundary, the global minimum turns into a maximum 
[Fig. [54|(a), curve (iv); Fig. [54|[c), curves (ii)-(iv)]. For the off-site modes, this maximum is always a global 
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(a) Off-site mode 


(c) 


On-site mode 



Figure 54: 

potential fi/U = 1.2 and the tunneling rates 2dJIU'. 
Ref. 13791 . ©2010 American Physical Society.) 


0.7 (i), 0.5 (ii), 0.3 (iii), 0.15 (iv), and 0.05 (v). (Reproduced from 


extremum. In the case of the on-site modes, the maximum of (h;) is either a global extremum [Fig. |54| [c), 
curve (iv)] or a local one which is accompanied by side minima [Fig. |54[c), curves (ii), (iii)[. Contrary 
to the results deep in the SF region, these types of the atomic distributions cannot be described by the 
DGPE. Similar features were also found for vortices |378L 142611427] and the underlying physical mechanism 
is essentially the same. 

The existence of the modes with maxima and minima of (n/) can be easier understood in the case of 
on-site modes. Let us assume that (h/) € (no — 0.5, ng -I- 0.5), where uq is a natural number including zero. 
Since = 0, the number of bosons at site Iq is fixed by some integer (local Mott insulator). From the 

minimization of the interaction energy, it follows that this integer coincides with ng which is either larger 
or smaller than (hj). This leads to the two types of solutions. For the off-site modes the situation is more 
involved but qualitative picture remains the same. 


In order to have a better understanding of the modes with the maxima of (h/), we depict in Fig. 55 


(/i, J)-diagram identifying the various types of solutions. The anomalous regions where (n/) attains a global 
maximum are almost the same for the off-site and on-site modes. They are, to a very good approximation. 


located in the “hole" areas of the (/r, J)-plane as displayed in Fig. 44 and hence, the corresponding modes 
can be interpreted as dark solitons of holes. The anomalous regions of the on-site modes which have minima 
of (ni) near the middle lattice point are located in the intermediate regions between particle and hole-areas 
and can thereby be interpreted as a mixture of dark solitons of holes and particles. With the increase of the 
filling factor the size of the MI lobes as well as of the anomalous regions decrease. 

In the MI phase, there is only trivial solution = 0, i.e., soliton modes do not exist. This follows from 
the fact that in the Gutzwiller ansatz, the excited states of the MI are products of local Fock states, where 
the occupation numbers n; can be locally different from the homogeneous filling ng. As a consequence, all 
must identically vanish and no soliton solutions are therefore possible within these parameter regimes. 

Experimentally, dark (or gray) solitons are typically created via a phase-imprinting method |42811429] . 
Initially (t = 0) the system of atoms is assumed to be in its ground state. During a short time timp one applies 
a spatially dependent potential on top of the lattice. In the Bose-Hubbard Hamiltonian, it is describ ed by 


the term time timp is much shorter than other characteristic time scales, from Eqs. (3231 

we get that the additional term induces a shift in the phase of the atomic states 


Qn(iimp) = exp(-i0/n) , 0i(fimp) = 0*-°^ exp (- 10 ;) 


(375) 


For the creation of dark solitons it is appropriate to choose a hyperbolic tangent imprinting potential, such 
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Onsite mode 



2dJ/U 


(b) 
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Figure 55: (a) Offsite modes. Dark (blue) areas show the first three MI zones (uq = 1,2,3). Gray (green) areas indicate 

the regi ons where the off-site modes have a global maximum of {hi) at the middle lattice sites around /q, see curve (iv) of 
Fig. |54| (a). In the rest part of the diagram, {hi) has only one extremum and takes the minimal value at the middle sites, see 
curves (i), (ii), (hi) of Fig. |54| (a). The lines of constant {h) corresponding to the ground-state densities are shown as well and 
labeled by the numerical values, (b) Onsite modes. Gray (green) areas depict the regions where the on-site modes have a 
global maximum of { hi) a t the middle lattice site /q, (in these regions {hi) have only one extremum which is a global maximum), 
see curve (iv) of Fig. |54| (c). In the light-gray (yellow) areas, the on-site modes have side minima near the maximum of {hi), 
see curves (ii) and (hi) of Fig. |54| (c). In the rest part of the diagram, {hi) has only one extremum and takes the minimal 
value at the middle site, see curve (i) of Fig. |54| (c). (Reproduced from Ref. [379], ©2010 American Physical Society.) (For 
interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 


that 


(t>l = 


e/ii, 


h 


A(j) 


1 + tanh 


l-lo 

0A51 


imp 


(376) 


where Iq is the middle point of the lattice. Here, /imp is the width of the interval around / = Iq where (pi/A(p 
grows from 0.1 to 0.9, and A(p is the amplitude of the imprinted phase [430]. Apart from the moving gray 
soliton, the phase imprinting also induces a density wave propagating in the opposite direction to the soliton, 
which appears due to the impulse imparted by the imprinting potential |428H43(J| . Numerical simulations 
performed in Ref. [575] according to this procedure show that the form of the propagating dark solitons 
becomes qualitatively different from those predicted by the DGPE if the system parameters are in the green 


(gray) regions of Fig. 55 where the standing solitons have global maxima of (h/). 


10.7.2. Stability of standing solitons 

We consider small perturbation of the soliton state determined by the coefficients c® as follows: Ci„(t) = 
where uji is given by Eq. (374) and 


cW(i) = u,„e— 


(377) 


Substituting this expression into the Gutzwiller equations and keeping only linear terms with respect to uin 
and vin, we obtain the system of linear equations: 


hujuin = 


E (^“7 


Ul'n' ^nl ^I'n' ) 5 


-hwVln = ^ (h”/ UVn’ + A”/ U//„/) 


(378) 
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Figure 56: Maximal imaginary part of the complex eigenvalues a; for the (a) off-site and (b) on-site modes. (Reproduced from 
Ref. |379| . ©2010 American Physical Society.) 


where 


An'l' 

^nl 


— n{n — 1) — iin — huji 


5l',lSn',n - J + 2(d - 5/'.; (v^ 5n',n+l + V^(5n,n' + l) 


- J 


Vn+1 Vn' + l c|°2+i 4°l'+i + Vn Vn' CjV;,_i 




[<5i,i' + l + (5z',/+l + 2(d — , 


TDn'l' 

^nl 


J Vn + l c|°2+i + VnVn' + l c|°Ai 4>'+i + <5i'.i+i + 2((i - . 


Eqs. (3781 generalize Eqs. (3381 to the inhomogeneous case and are analogous to the Bogoliubov-de Gennes 
equations which were employed for the stability analysis of the dark solitons governed by the DGPE |3811 

The stationary modes are linearly stable, if all the eigenvalues hu) are real. Numerical solutions of 


the eigenvalue problem (3781 show that most of the eigenvalues are real but there are always few ones, 


which contain nonvanishing imaginary part uii- The magnitude of uii determines the inverse lifetime of the 
solitons, which can be almost equal or drastically different for the off-site and on-site modes and there is no 
any principal difference in this respect between the normal and anomalous modes. 

Figure shows the maximal imaginary part of the complex eigenvalues w which vanishes in the MI 
regions, where solitons do not exist, but does not vanish in the SF region. With the increase of /i and J, 
maximal uji increases for both types of soliton modes meaning that the instability grows. There is, however, 
one important qualitative difference between the off-site and on-site modes. For the on-site modes, there 
are rather small regions between the MI lobes, where uJi is close to zero and much smaller than that for 
the off-site modes, i.e., the on-site solitons are much more stable. This feature has some similarity to the 
stability of the standing dark solitons governed by the DGPE, where it was found |381| that on-site modes 
are stable if the tunneling J does not exceed a certain critical value, while off-site modes are unstable for 
all tunnelings. 


11. Spinless bosons near Feshbach resonance 

Magnetic Feshbach resonances provide and important tool to tune the interaction strength in ultracold 
atomic gases. They found numerous applications and are often used in experiments for the observation 
of different phenomena, some of them were already mentioned in this review. Feshbach resonances occur 
when the energy of two scattering atoms is close to the energy of their molecular bound state. In essence, 
it is a quantum interference effect of these two states. The energy difference between the states can be 
controlled by an external magnetic field when the corresponding magnetic moments are different which 
gives a possibility of tuning. Excellent reviews on Feshbach resonances in the context of ultracold atoms 
were given in Refs. |34ljMl l43TII433| . However, the lattice problems were not addressed there and the aim 
of the present section is to fill this gap. 
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11.1. Hamiltonian 

We consider a system of identical atoms, each having the mass M, in an optical lattice created by a 
far-detuned standing laser wave. In addition, the atoms are subject to the external magnetic field B, with 
B = Bq corresponding to the Feshbach resonance of the width AB. This allows to convert pairs of atoms 
into molecule of mass 2M and there is a reverse process when a molecule dissociates into two atoms. In 
order to describe Feshbach resonance we adopt two-channel model with one open channel that describes 
asymptotically free atoms and one closed channel corresponding to the two-atom molecule [Ml US]. The 
model is adequate to describe an isolated (narrow) resonance. 

The far-detuned optical laser field creates an effective potential not only for the atoms but also for the 
molecules that has the same form as for the atoms but the amplitude is doubled. Using the basis of the 
atomic and molecular Wannier functions Wa(x), Wm(x) one can derive the discrete Hamiltonian of the 
system. In the tight-binding approximation, the microscopic Hamiltonian is given by |485H487] 


- Jy 


a 

1 I/ = l 


h.C.) - 


1 U=1 




u 


afalojOj -I- g 


, (379) 

1 


where dj (6|) and dj (Sj) are creation and annihilation operators of a single atom (molecule) at a lattice 
site 1, d = Ag,{B — Bq) is a detuning from the Feshbach resonance. Here, Ag, is the difference in magnetic 
moments of the two atoms and a molecule. The atom-molecule conversion is determined by USD 


g = H 


2TrasABAii 


M 


IUf(x)IFm(x)d> 


(380) 


where IFa(x) and IFni(x) are the Wannier functions for the atoms and molecules, respectively. In the 
three-dimensional lattice and in the Gaussian approximation [see Eq. @1, it has the form [43511438| 


g = h 


2'KasABAg, 


(381) 


The on-site molecule-molecule and the atom-molecule interaction parameters C/m, Uam are defined by the 
expressions similar to that for the atomic parameter U. Due to the differences in the physical properties of 
the atoms and molecules discussed above, the molecular tunneling parameter is much smaller than the 


atomic one. The Hamiltonian (3791 does not conserve the number of atoms and molecules separately but 


the total number of the atomic constituents 


iVt 


=i:».=i:( 


Oj dj 


26j'6i 


(382) 


is preserved. 

In the regime of small tunneling and in the case of two atomic constituents, the low-energy properties of 
the system can be described restricting the local Hilbert space of the Hamiltonian (3791 by the states with 
two atoms or one molecule. This leads to the mapping on the spin-1/2 quantum Ising model [439] : 


a 

-^ani ~ Rising — Jz EEsfs.+..+E( 


1 i/=i 


const . 


(383) 


The spin operators are given by 

waa 


S+ = 


V2 


S- = 


a^a^b 

~VT 


5" = 
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s+ + s- 
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y^b — d^d/2 


(384) 
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Figure 57: Eigenenergies in the case of two atoms on the same lattice site. Solid lines show the results given by Eq. ( |387[ l 
which correspond to the lowest-band approximation. The results for the infinite number of bands [Eq. ( |390[ l] are shown by 
dashed lines. U = 0, 2y/7rg‘^l = O-l- 
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where we omitted the site index 1. The first term in the Hamiltonian (3831 stems from the virtual hoppings 
of atoms and molecules between the neighboring sites and corresponds to an effective magnetic exchange 
interaction with the coupling constant 


Jz = 


4J2 


- U 




Other parameters 


hii=S-U, hj_ = 2gy/2 


(385) 


(386) 


correspond to the longitudinal and transverse magnetic fields. This mapping facilitates understanding of 
the phase diagram of the system which will be discussed in Section |11.4| 

11.2. Two atomic constituents on one lattice site 


In the limit of vanishing tunneling, the Hamiltonian (3791 becomes a sum of local terms. Then the on-site 


problem for two atomic constituents can be easily solved analytically. In this case there are two eigenmodes 
which are superpositions of the two-atom and molecular states with the energies 




and the probability to find a molecule is given by 


'S-U 


+ 2^2 , 


(387) 


Pm± - 2 


1± 


S-U 


^{5-Uf + S~g\ 


(388) 


If the system is initially prepared in the state with two atoms at each site, the probability to find a 
molecule will oscillate in time according to 


4^^ 

Pm± = ^ (1 - COSCji) , 


Ej- — E— 


^{5-Uf + ^g^ 


(389) 


This kind of Rabi oscillations was observed in the experiments with ®’^Rb in deep optical lattices |438| . 

It is interesting to compare this solution with that obtained in Ref. |435| for two atomic constituents 
on one site of a deep lattice for the infinite number of bands neglecting the atom-atom interaction. In this 
case, each site can be described by a harmonic potential with the frequency Swho and the eigenenergies E 
are determined by the equation 


^ V{-E/2njw^,) 

/iWho T[-E/2Eo-i,o - 1/2) 
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where r(a;) is the Gamma-function. The eigenenergies given by Eq. (387) for [7 = 0 and the solutions of 


Eq. (390) are plotted in Fig. 57 As we see, the lower branch in Eq. (387) is in excellent agreement with 


the lowest branch of Eq. (390) for arbitrary S. The upper branch fails to reproduce the second branch 


of Eq. (390) if 5 is far above the Feshbach resonance where the contribution of the second band becomes 
significant, remaining however in a very good agreement near the resonance and below it. This implies 
that the lowest-band approximation is valid if the detuning S is less than the gap between the two lowest 
Bloch bands, which is the quantity of the order of hoJho, and/or if we are interested in the eigenmodes of 
the Hamiltonian (379) with the energies less than the energy of the second Bloch band. In addition, the 


parameters U and g must be much smaller than the bands separation. 

11.3. Two-body eigenmodes and bound states 

As in section [63| we consider a one-dimensional model with L lattice sites assuming that L is odd. Under 


periodic boundary conditions the eigenstates of the Hamiltonian (379) for two atomic constituents have the 
form 

= (391) 


7=0 


tk 


L-l 


where | ATI) is the atomic part defined by Eqs. (187), (188) and |li„0 ... 0) is a state with one molecule on the 


first lattice site and all the other sites being unoccupied. The eigenvalue problem for the Hamiltonian (379) 
can be written down as follows: 


^kc^q, + V^gcKCio — , 

T ckuq + H^ckq.1 = E^^oxna ) 

(L-l)/2 

CKnr' = E^^ckut , T = 1, • • ■ > 

r'=o 


(392) 


L-l 


where 5k = 6- 2 Jya cos {Ka) and the nonvanishing matrix elements Hk^ are determined by Eq. (^190|. 
The normalization condition takes the form 




(L-l)/2 

E 

r=o 


Ici^nr 


= 1 . 


(393) 


The aim of the present section is to investigate the influence of the molecular mode on the bound state |4401 - 

\m\ . 


As in the case of two atoms considered in section 6.3 the spectrum consists of two types of modes. The 


scattering modes have the same energies as in Eq. (153). The bound states have the form (195). Substituting 


this ansatz into Eq. (392), in the limit L —>■ oo we obtain the equation for the eigenenergy 

2^2 


Sr — Ur E Qk , Ur — U -\- 


£k — 5r 


(394) 


where gR is defined in Eq. (153). Formal comparison with Eq. (196) shows that Ur plays a role of the 


effective atomic interaction. The corresponding values of crq and bR are given by 

Ur — Sr 


cro = 


{£r - SrY (1 - bl) 


{£k - SrY (1 + bl) + 2g2 (1 _ bl) 

and the amplitude of the molecular state takes the form 

V^gcRo 
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QK 


£r - Sr 
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Figure 58: Energy spectrum of two atoms in the presence of Feshbach resonance. IJjJ = —5, (5/J = —4.5; gjJ = 0 (a), 
S:! 1.118 (b), 2 (c), S:! 3.35 (d), 4.15 (e), 6 (f). 


Eqs. (394), (395) remain unchanged under transformation U —i —U, Sk —t —(5ic, £k —£kj bx —bx- 


Therefore, in order to get the complete solution it is enough to study the case of attractive interaction in 
the whole range of Sx and g. The latter can be considered as positive because Eq. (394) contains only. 


Eq. (394) can be multiplied by {£x — Sx) and treated as quartic equation for £x which contains always 


four roots. However, depending on the values of the parameters only one or two roots are real and provide 
normalized eigenstates with \bx\ < 1 implying that the others are unphysical and should be rejected. In the 


special case g = 0, vanishes and Eqs. (394), (395) reduce to the solution (196), (195) in the absence of 
the Feshbach resonance. 

Although analytical solutions of the quartic equation are well known, simple expressions for £x can be 
obtained only in some special cases. For instance, one can easily show that in the special case Ja = Jm = 0 


the physical solutions of Eq. (394) are given by (387). 


In the limit of large detuning, |5ic| 3> \£k\-, the effective interaction parameter takes the form 

Ux = U — 2g‘^jSx 

which is equivalent to the expression for the effective scattering length 

= as(l — ABAfi/S) 


(397) 


(398) 


that appears in the mean-field theory as a result of the adiabatic elimination of the molecular field |445| . 
In this limit, the solution for £x is given by Eq. (196) with U being replaced by Ux- In the special case 
(7 = i5ic = 0, we obtain 


£x± = ± 



(399) 


In general, Eq. (394) allows at least one bound state which corresponds to that for the two atoms without 


Feshbach resonance in the limit of vanishing g. With the increase of g, the absolute value of the energy of 
this state grows but the sign remains unchanged. The second bound state exists, provided that the absolute 
value of the quantity Qx = {USx — 2^^) / (Uqx) is larger than one. Qx > 1 leads to positive Ux and 
Qx < —1 gives negative Ux- As a consequence, the second mode exists for any value of K, if jQol > 1- 
Otherwise, the second mode exists only for \K\ > Kc = ^ arccos jQol within the first Brillouin zone. In the 
expression for we have neglected the molecular tunneling J^. For Qq = 0, = tt/o, i.e., the second 

bound state does not exist. Different types of solutions are shown in Fig. 

The properties of the bound states are determined by the values of the effective interaction parameter. 
The state is attractively bound for Ux < 0, even if the interaction parameter U is positive or vanishes. If 
Ux > 0, the state is repulsively bound, no matter what the parameter U is. Moreover, attractively and 
repulsively bound states can coexist which is not possible in the absence of Feshbach resonance. 
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Figure 59: Phase diagram in the local limit (Ja = Jm = 0) vs chemical potential fi and detuning (5 from the Feshbach 
resonance, g = —0.5, U = Um = !• Uam = 1 (left), Uam = 0.25 (right). Each phase is labeled by the number no of atomic 
constituents per site. (Adapted with permission from Ref. l4i6l V 


11.4-- Zero-temperature phase diagram 

In the absence of tunnelings, Ja = Jm = 0, and for an arbitrary integer number (n,) = tiq of the atomic 
constituents at each lattice site, the local eigenstates of the Hamiltonian are superpositions of the form 

[no/2] 

1//’)= X! C'„„|na,nm) , (400) 

n-m—O 


where is the number of molecules and the number of atoms = uq — 2n^. The (/i, 6) - diagram worked 
out in Ref. |446| by numerical diagonalization of the Hamiltonian matrix for different values of uq is shown 
in Fig. 59 It consists of MI regions labeled by the values of hq. If the detuning is large {\5/g\ 1) and no is 

even, the MI state exists on both sides of the resonance. For positive <5 it is dominated by \n^ = no, nm = 0) 
and for negative 6 by \n^ = 0,nm = no/2). For odd no and still for large \&/g\, the state exists only 
for positive <5. For large negative 5, the ground state for odd no is a superposition of Fock states with 
Ura = {no ± l)/2 and is expected to be unstable against superfluidity for any finite tunneling. The phase 
diagram for small detuning strongly depends on Uam'. For small C/am, the MI phase with odd no can appear 
also for negative 5. 

In order to understand the phase diagram for nonvanishing tunneling parameters, it is useful to pay 


attention to the fact that the Hamiltonian (3791 is invariant under U(l)x 2'2 transformation 1 44711449] : 


h\ —>■ h\e 


i29 


ai 


±aie 


ie 


(401) 


where 9 is real. In the dimensions larger than one, the continuous U(l) symmetry may be broken without 
breaking the discrete Z 2 symmetry Si —±Si. This leads to the molecular condensate phase (MC) with the 
order parameters ( 61 ) yf 0, (Si) = 0 and corresponds to the Ising degree of freedom in the disordered phase 
coexisting with the molecular superfluidity. On the other hand, the U(l)xZ 2 symmetry can be completely 
broken which gives rise to the atomic-molecular condensate (AC+MC) with nonvanishing (Si), (Si) and 
corresponds to a Z 2 ordered Ising degree of freedom coexisting with the atomic and molecular superfluidity. 
If the U(l) symmetry is not broken, the system will be an insulator as in the case of one-component bosons. 
Mean-field studies reported in Refs. |4351144611450] for commensurate and incommensurate fillings are in 
agreement with these general considerations. 

In one dimension, breaking of the U(l) symmetry is prohibited and the formation of nonvanishing 
expectation values (oi) and {b\) is excluded. Low-energy effective theory based on the “bosonization" ap¬ 
proach |281j shows that in the atomic-molecular superfluid phase which is analogous to AC+MC in higher 
dimensions the atomic and molecular one-body correlation functions have a power-law decay at large dis¬ 
tances mzHiis]: 

{ala,J ^ |/i - /2|-“- , (bib,^) - 1^1 - . (402) 


Due to the phase locking of the atomic and molecular components arising from the Feshbach term in the 
Hamiltonian (3791, = dofa- Fr the other superfluid phase corresponding to MC, the molecular correlation 
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Figure 60: (color online) Phase diagram of the ID Hamiltonian ( |379[ | with two atomic constituents per site, showing a Mott 
insulator (MI), a molecular condensate (MC), and a coupled atomic-molecular condensate (AC+MC). The parameters are 
Jjjj = Ja/2, g = t//2, Um. = Uam = U/A. The squares and circles indicate the vanishing of the one-particle and two-particle 
gaps, Eig and E 2 g, respectively, as L ^ oo. The stars and crosses indicate where the molecular and atomic correlation 
exponents, Om and reach 1/4 in the MC and AC+MC phases respectively. These values correspond to a molecular and an 
atomic Berezinskii-Kosterlitz-Thouless transition, respectively. (Adapted with permission from Ref. | 449 |. ©2012, American 
Physical Society.) 


function has the same power-law decay but the atomic function decays exponentially |448l I449| 




(403) 


where ^ is the Ising correlation length. However, the correlation function of atomic pairs exhibits a power-law 
behavior |448L 1449] : 

~ 1^1 “ 4r“” ■ (404) 

In both superfluid phases, the atomic-molecular correlation function decays as a power law with 

the exponent and all particle-number correlations have the same asymptotics 


(^a^i^/3^2) ~ (^a^i ) (^/3£2 ) 


Cap 


(405) 


where a, /? = a, m and Cap are nonuniversal constants. 

These field-theory predictions were successfully tested by DMRG calculations in one-dimensional chains 
up to L = 512 sites with two atomic constituents per site |448L 1449] . It was found that MC (AC+MC) 
undergoes MI transition, when the molecular (atomic) correlation exponent am (aa) reaches the value 1/4 
which is exactly the critical value for the Berezinskii-Kosterlitz-Thouless transition. At the transition point 
from the AC+MC phase into the MI, the molecular exponent takes the value am = 1- This is an indication 
of the molecular superfluidity which signals the absence of a single-component atomic superfluidity close to 
the MI boundary in contrast to Ref. |451| which was claiming the opposite based on QMC calculations. The 
latter was attributed to finite-size effects mg. 

The transition between the two SF phases is expected to be in the universality class of the (d + 1)- 
dimensional Ising model. The values of the critical exponents of the correlation length and of the order 
parameter for the two-dimensional Ising model {v = 1, f3 = 1/8) are in perfect agreement with the DMRG 
calculations in one dimension |44811449| . 

The phase diagram of the one-dimensional chain with two atomic constituents per site obtained in 


Ref. |448| by DMRG method (see also Ref. |449j l is shown in Fig. 60 The phase boundaries correspond 
to the vanishing of one-particle and two-particle excitation gaps (n = 1,2) Eng = p,n+{Nt) — /i„_(iVt) 
extrapolated to the thermodynamic limit, where 


/i„±(iVt) = ± [EoiN, ± n) - Eo{Nt)] jn 
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and i?o(-/Vt) is the ground-state energy for a system of fixed size L with the total number Nt of the atomic 
constituents. In the second-order of strong-coupling expansion, the MI phase reveals a phase transition from 
Z 2 disordered phase with vanishing staggered magnetization {Sf)/L, where = (hra — ns.l‘Z) 12, to 

the ordered phase with a finite staggered magnetization and long-range antiferromagnetic correlations |439| . 

The phase diagram similar to that shown in Fig. 60 was also obtained in Ref. m by means of QMC 
simulations. However, the MI phase was named “super Mott" (see also |452p based on the observation that 
the superfiuid stiffnesses of the atomic and molecular components calculated from the fluctuations of the 
corresponding winding numbers do not vanish, although the superfiuid stiffness of the whole system does. 
This interpretation was criticized as having no sense due to violation of the particle-number conservation in 
the atomic and molecular subsystems |453| . 


Zero-temperature phase diagram of the Hamiltonian (3791 was recently studied in the case of two atomic 


constituents per site for two-dimensional and three-dimensional lattices within the framework of the de¬ 


coupling mean-field approximation [see Eq. (318l] complemented by QMC calculations for two-dimensional 
lattices |450| . These studies reveal the same superfiuid phases as in one dimension accompanied by true 
Bose-Einstein condensation in agreement with the previous works |4351 \4A6\ . It was pointed out that one 
has to distinguish between three types of insulators: molecular, atomic-molecular and Feshbach insulator. 
While the molecular insulator is a conventional one with the same properties as in the case of one-component 
spinless bosons, the atomic-molecular and Feshbach insulators are basically products of the superpositions 
of local states with two atoms and one molecule on each site. The main difference between the atomic- 
molecular and Feshbach insulators is that the two-particle excitation gap E 2 g decreases with the increase 
of g for the former but increases for the latter starting from E 2 g = 0 at 5 = 0 and the transition between 
the two insulators is rather a crossover. The transition from the superfiuid into the Feshbach insulator as 
well as the transitions from the atomic-molecular superfiuid to all the insulating phases were found to be 
strongly first order. It was also established that the mean-field theory captures correctly the succession of 
phases in the system. 

In the case of two-component bosons corresponding to different atomic species near the Feshbach res¬ 
onance which supports the creation of heteronuclear molecules, the zero-temperature phase diagram has a 
reacher structure |454L 1455] . It contains MI regions surrounded by the regions of single-component atomic 
and molecular superfiuids (in contrast to the homonuclear case discussed above) as well as the regions with 
all three superfiuids. However, there are no phases with just two superfiuid components. In the regime 
of small tunneling with two atomic constituents per site, the underlying Hamiltonian can be mapped to 
the quantum Ising model with longitudinal and transverse fields as in the homonuclear case and, therefore, 
reveals again the Ising phase transition. 

Up to now the complete phase diagram of ultracold atoms in lattices near Feshbach resonance was not 
explored experimentally. So far only the creation of the homonuclear molecular Mott insulator has been 
reported |45611457] . 


12. Spin-1 bosons 

If the atoms a trapped by purely optical means, the spin degree of freedom is not frozen. We consider a 
dilute gas of bosonic atoms with hyperfine spin F = 1 possessing three Zeeman-degenerate internal ground 
states with magnetic quantum numbers a = mp = 0, ±1 in the field of an optical laser described by a 3 x 3 
matrix U^®'®(x). The system is governed by the following Hamiltonian |458l 1459] : 


Hf=i = 




(407) 


+ ^^l(x)^j3(x)4'/j(x)^„(x) -k ^4'),(x)^|^,(x)F„^ • F„/^/^;3 /(x)4';3(x) 
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where 'l'ct(x) is the bosonic field annihilation operator for the atom in the hyperfine ground state j-F = 1, a). 
F is a vector of traceless spin-1 matrices: 


0 1 0 
= 1 0 1 
v2 \ 0 1 0 




0-10 
0 -1 

1 0 



(408) 


In Eq. (4071, summation over repeated spin indices a,/3 is implied. 


The parameters gs,a describe the strength of the repulsive interactions of the atoms and the spin-changing 
collisions. In three dimensions, they are related to the scattering lengths oq and 02 of two colliding bosons 
of mass M with total angular momenta 0 and 2 (singlet and quintuplet channels) as 


9s = ^ (ao + 202 ) , g,= —{a 2 - ao) . 


3M 


3M 


(409) 


The values of og and 02 for atoms usually used in the experiments have been reviewed in Refs. |46011461] . 
In the case of ^Li, and ®^Rb, g^, turns out to be negative, while for ^^Na it is positive. 

As long as all the atoms are in the internal state with a = —1 or a = -1-1, the spin degrees of freedom 
do not play any role and all the physical properties will be the same as for spinless atoms. The spin degrees 
of freedom come into play when the internal states with different values of a are populated by many atoms. 
This property can be used to probe the particle-number statistics across the MI-SF transition |462| . If the 
atoms are initially prepared in the internal state \F = l,mF = —1) and then transferred to the state with 
mp = 0, a, reversible exchange between the populations in the mp = 0 and mp = ±1 Zeeman sublevels is 
observed unless the initial state is a product of local Fock states with the occupation numbers equal to one. 


12.1. Bose-Hubbard model 

If the detuning of all the lasers creating a n op tical lattice is much larger than the fine splitting of the 
electronic energy levels, matrix F'®'® in Eq. (4071 becomes a scalar, i.e., all spin-components a = 0,±1 


experience the same lattice potential as in the case of spinless atoms. Expanding the field operators in the 
Wannier basis and using the tight-binding approximation, we obtain |9I[ 146011463H465] 


= -J 


a 

EE( 


«ia«l+e„.a+h-C. 


1 




yhi (ni - 1) -k Y (tf - 2hi 


(410) 


where 

is an operator of the total number of atoms on site 1 and 

Li = aj^Fa/^di^ 


(411) 

(412) 


is the spin operator on site 1. Its components obey the standard commutation relations for the angular 
momentum 

L\a\ ; F\a2 


= lea 


3 Lias 


(413) 


where e is the completely antisymmetric Levi-Civita tensor. T he op erator fii commutes with Li, and the 
total spin operator L = Li commutes with the Hamiltonian (4101. 

The interaction parameters Us a are given by Eq. (551 with g replaced by gs a and their ratio is limited 

by 

-I<|<i. (414) 
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provided that both ag and 02 are positive. The tunneling matrix element for the nearest-neighboring sites J 
is exactly the same as in the case of spinless atoms. The term of the Hamiltonian (4101 with the prefactor 
t/a has an explicit form (the site index is omitted) 


- 2 h = 


a\a\diai + ~ 2d\a}_idid_i + 2a|oJa]^ag 


-I- 2d}_^dl,d_idQ -I- 2agaja]^a_]^ -|- 2a^2^a|aQag , 


(415) 


where the last two summands describe the spin-changing collisions. The latters were observed in experiments 
with ®^Rb atoms in deep optical lattices |466L 1467] and the measured differences of the scattering lengths 
qq — 02 agree with the theoretical predictions |468| within 20%. 

Sometimes it is convenient to work with the operators ba, a = 1, 2, 3, |465L 1469] 


h = 


a_i — ai 


V2 


62 = - 


.a_i+ di 
^— 

72 


bs — do , 


(416) 


which satisfy the standard bosonic commutation relations and transform as vectors under spin rotations. In 
therms of the operators ba, the particle-number operator can be expressed in a usual way, n\ = bla^a’ 
the spin operator is given by 


L 


lap 


— *^010203^1026103 ) 

= hi (hi -I- 1) - b[^M. b,,._b. 


^lai'^lai *^102 ^ 1^2 * (417) 

The transformation properties of ba, a = 1,2, 3, can be verified with the aid of the commutation relations 


La\ , ba 


= le. 


010203^0,3 5 


^01 ; 6^2 


— *£010203603 . 


Using (4161 and (4171 the Hamiltonian (4101 can be rewritten in the form 

d 


z/=l 1 


^ Lw . o + h - c .) + Y. 
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U, + Ua 


Ua 


hi {ni - 1) - ^blMa.ha^ha^ 


(418) 


(419) 


which is invariant under global spin rotations. 

As in the case of spinless bosons, the eigenstates of the spin-1 lattice system can be studied using the 
basis of local Fock states jnn, nio, ni_i). On the other hand, magnetic properties are better described in the 
basis of spin states |ni, Li, L 13 ) where the quantum numbers label the eigenstates of the operators hi and Li. 
The spin states can be uniquely expressed in terms of the Fock states |47()11471] and the Hamiltonian (4101 
imposes two constraints on the quantum number Ly. First, the total spin cannot be larger than the total 
number of particles, i.e., Li < n\. Second, bosonic symmetry under permutation of any two particles leads 
to the requirement that ni and Li should have the same parity, i.e., Li must be even (odd), if ni is even 
(odd). A rigorous proof of this statement was given in Ref. |470| . 


12.2. Single-particle states 

Single-particle eigenstates in a homogeneous lattice with periodic boundary conditions can be written 
as |ka) = hj(.^|0), a = 0, ±1. Although being a trivial generalization of the spinless case, they show 
already different magnetic properties. |k0) is a simplest example of a nematic state that has vanishing 
expectation values of all spin-components, i.e., (L„) = 0 for a = 1,2,3, but breaks the spin symmetry 
because (L?) = = 1 and (L§) = 0 . 
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12.3. Eigenstates of two atoms 

We consider first the eigenstates of two atoms on a single lattice site. In this case there are six states 
and four of them are not influenced by the spin-changing collisions. In terms of the Fock states the latters 
are given by 


|V'i) = |2,0,0), IV'a) = 10,0,2), |V^3) = | 1,1,0), IV' 4 ) = |0,1,1) 

(420) 

Two other eigenstates have the form 


IV's) = y||0,2,0) + y||l,0,l) 

(421) 

IV'e) = y||0,2,0)-y||l,0,l). 

(422) 

In terms of the spin states, the six eigenstates are given by 


IV'i) = |2, 2,2) , 1 ^ 2 ) = |2, 2, -2) , IV- 3 ) = |2,2,1) , 

IV'4) = |2, 2, -1) , IV's) = |2,2,0) , IV-e) = |2,0,0) . 

(423) 

The state j'l/'e) is unique and has the energy E = Us — 2t4, while all the others are degenerate and their 
energy is E = Us-\- U^,. For positive C4, \i/g) is the ground state and for Hs, < 0 it becomes an excited state. 

Note that Itfe) is the only state among the others which has equal populations of all spin components 
a = 0, ±1, i.e., (hi) = (ho) = (n_i) = 2/3. It is a spin singlet and Eq. (4221 defines the creation operator 
of a sinelet oair |4f)5l I472| 

^Ig = (44 - 2alaLi) = ■ 

(424) 


Since commutes with the spin operator L 
a = 1,2,3. Moreover 


it changes neither the total spin nor the spin components 
from Eq. (4171 one can see that the eigenstates of the total spin operator must be 
always eigenstates of the “singlet counting operator" 

The eigenstates of two atoms in the case of nonvanishing tunneling can be readily constructed from the 
solutions for distinguishable and indistinguishable atoms discussed in section |6.3| In the present case there 
are six bound states with the effective interaction parameters U = Us + 11^, and U = Us — 2C4,. Due to the 
restriction dini , U is always positive and the energies of all bound states appear to be above the scattering 
(quasi-)continuum. 


12 . 4 . Ground-state phase diagram 

In this section we shall discuss the ground-state phase diagram of the lattice spin-1 system in the absence 
of an external magnetic held. As in the case of spinless bosons, it consists of the MI and SF phases. However, 
the spin degree of freedom leads to new interesting aspects related to magnetic properties which appear to 
be completely different in the case of positive and negative Ug_. 

The difference between positive and negative 14, can be already expected looking at the ground state of 
the Hamiltonian (4101 in the limit of small tunneling. In the case of negative U^ the ground state should 
prefer the largest possible values of Li, while in the case of positive t/a the smallest possible values of L\ 
are favorable. Moreover, we can expect differences between even and odd fillings in the case of positive 14 
because the smallest value of L\ is zero in the former case and one in the latter. 


12.4.1. Ha = 0 

In the special case L4 = 0, the spin-dependent interaction is absent and the numbers of bosons in all 
components a = 0, ±1 are conserved separately. The ground states are highly degenerate and exhibit “SU(3)- 
ferromagnetism " ms]. In the case of integer fillings and for small J/Us, the ground state is predicted to 
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be a nematic insulator |463| . This state has vanishing expectation values of all spin-components: (La) = 0, 
a = 1,2,3, but the spin-rotational symmetry SO(3) is broken, while the time reversal symmetry is preserved. 
The order parameter for the nematic state is a traceless symmetric tensor Q with the entries |4631146511472] 

Qab = (£aib)-^(L2), a, 6 = 1,2,3. (425) 

In the case of integer fillings and large J/Us as well as in the case of fractional fillings and arbitrary 
J/Us, the ground state is a polar superfluid. In two and three dimensions, it breaks both U(l) and SO(3) 
symmetries ssa. 

Mean-field theory predicts the same zero-temperature phase diagram as in the spinless case and, in 
particular, continuous transition from the MI into SF |474] which agrees with QMC calculations in two 
dimension mg. However, at small but finite temperatures, the transition becomes first order and if the 
temperature is increased further it is again continuous |474] . 


12.4.2. t/a < 0 

According to the theorem proven in Ref. \m\ . the ground state in the case of negative Ug, exhibits 
saturated ferromagnetism {L\ takes the largest possible value) in any dimension and both in the MI and SF 
phases. Mean-field theory predicts in this case continuous transitions form the MI into SF with the phase 
boundary described by Eq. ( |329 1, where U = Us + Ug, |476| . Continuous character of the phase transition 
was confirmed by QMC calculations in one |477| and two dimensions |475| . It was also shown that the phase 
diagram in one dimension can be described very accurately using the third-order strong-coupling expansion 
for spinless bosons with U = Us + U^ |477| . In two dimensions, the MI lobe for (ni) = 2 obtained by the QMC 
method is well reproduced by the mean-field result, whereas the MI lobe for {n\) = 1 in QMC calculations 
is significantly larger than in the mean-field theory |475| . QMC calculations also demonstrated that the 
populations of the spin-components in the MI and SF phases are (hio) = 2(nii) = 2(hi_i) = (ni)/2 |475| in 
agreement with the mean-field theory |474| . Note that in the case of two atoms on one lattice site the state 
\ip 5 ) determined by Eqs. (4211, (4231, which has the largest possible spin L = 2 and vanishing L3, shows the 
same ratio of populations. 


12.4.3. [/a > 0 

In the case of positive f7a th® system has a rich variety of phases with different magnetic ordering and 
the dimensionality plays an important role. We consider first the case of even integer filling (hi). In the 
limit of vanishing tunneling, the ground state for each isolated site is a spin singlet (which is unique) and 
the excitations are spin Li = 2,4,..., (hi) which are gapped by the energies of the order of [4,. Therefore, 
in all dimensions the ground state in this limit is a spin singlet Mott insulator |478| . QMC calculations in 
one and two dimensions for (hi) = 2 show that the populations of the components a = 0,±1 are equal to 
each other within this phase |47511477] which is a general property of the singlet states |461| . 

If the tunneling grows, the spin singlet state transforms into the insulating nematic state |463L I465L14721 
I475L1477] . provided that the ratio U^/Us is small enough |465L 1475] . Mean-field theory predicts in this case a 
first-order phase transition |46511472| and provides an estimate of the transition point Jc which depends on 
the filling. For two particles per site, J/ = f7st7a/(4d), and for large even fillings j"/ = 9C/sf7a/(2d(hi)^) [465] . 
The insulating nematic state exists only for small enough values of U^/Ug, which are less than 0.05/d in 
the case of two particles per site. These predictions of the mean-field theory were confirmed by QMC 
calculations in two dimensions m- However, in one dimension QMC calculations reveal that the advent 
of the insulating nematic state is a crossover [TTZ] . 

If the tunneling is increased further, the system undergoes a transition into the SF phase. However, the 
character of the transition depends on Ug^/Us'. if it is small, the transition is first order but for larger values it 
becomes second order. This result was obtained first within the framework of the mean-field theory |4741I479] 
and then confirmed by QMC calculations in one and two dimensions |475L 1477] . Earlier mean-field studies 
provided an estimate of the largest value of Ua,/Us which still allows the first-order transition in the case 
of (fii) = 2: Ug,/Us 0.32 |479| . which was later corrected to C/a/tA 0.2 |476| . The latter is in a better 
agreement with QMC calculations in two dimensions: Ug,/Us ~ 0.15 |475| . 
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If {n\) is odd, the situation becomes different. In the absence of tunneling, the ground state for each 
individual site is a spin L\ = 1 state with a three-fold degeneracy. For two decoupled sites, the ground 
state has a nine-fold degeneracy corresponding to the states with the total spin Ltot = 0,1,2. Finite 
tunneling lifts this degeneracy but the form of the resulting state depends on the dimensionality and other 

details IMIIMIIIZH]. 

In one dimension, the MI phase in the case of odd fillings is always dimerized |46511478L148014485] . The 
dimer state breaks translational symmetry and favors singlets on every second bond. Its distinguishing 
property is the doubling of the unit cell of the lattice, while the spin long-range order is absent. The 
simplest dimerized state can be written as |465L 1482] 

\D) = (g) = 1, L,+i = 1, + L^+i = 0) , (426) 

i odd 

where the product could be also over even i, i.e., the state is doubly degenerate. The dimerization can 
be described by looking at the expectation values of a pair Hamiltonian adjacent bonds {H = 

The corresponding order parameter reads [48111482] . 

If the parameter [/a/C4 tends to zero but remains finite, the amplitude of the dimer state becomes 
very small |485| and the spectrum of excitations shows qualitative changes |48611487] . This might be an 
indication that the MI state in a one-dimensional lattice becomes nematic. However, the limitations of 
different analytical and numerical methods applied in this regime do not allow to make a certain statement 
(see, e.g., discussion in Ref. |485p . In higher dimensions, insulating phases with an odd number of particles 
per site are always nematic |46511472[ 14751 IT55] . 

With the increase of the tunneling parameter, the system enters again into the SF regime. For (hi) = 1, 
the transition is always continuous |47444477l 1479] . For larger odd fillings and small values of Ua/Us, it can 
also become first order, although in this case the effect is not so pronounced as for even fillings gzni. 

The formation of the singlet pairs in the case of even fillings stabilizes the MI phase against the transition 
into the SF. This leads to the asymmetry of the phase diagram: the size of the MI lobes is larger than for odd 
fillings. This feature was demonstrated by DMRG |48144484] and QMC calculations |475] and also captured 
by the mean-field theory |4711147611479] . 


12.5. Effective spin-ll2 Bose-Hubbard model 

Lin-0-lin laser configuration discussed in section 12.4] leads to two sets of orthogonal Bloch eigenmodes 
denoted by the indices 0 and A. This allows to derive effective spin-1/2 Bose-Hubbard model from the 
Hamiltonian (4071. In the tight-binding regime, the atoms stay always in the lowest Bloch bands with the 


dispersion relations ifQ°^(k) and £’g^'^'’(k). Then the spinor-field operator '5'(x) can be decomposed as 


AA). 


’5'(x) = 


1 cr=0,A 


Wj")(x)a,,i, 


(427) 


where Oo-i is the Bose annihilation operator for the cr-mode attached to the 1th lattice site. w|‘^^(x) = 
W(‘")(x — xi) are three-component Wannier spinors for the lowest energy bands localized at the minima 
of the lattice potential labeled by 1, which have the form similar to Eq. ( ]40| . They are obtained by the 
solution of the eigenvalue problem for the single atom discussed in section 12.4] and satisfy the orthonormality 
condition 

J W|f^^(x) • )(x) dx = . (428) 


Substituting Eq. (4271 into Eq. (4071 and taking into account only the hopping between the nearest lattice 


sites and the on-site atomic interactions, we obtain the two-component Bose-Hubbard Hamiltonian |126] 


Hbh = (®ai«<Ti+e,. + h-c.) - 1) 


+^E 


noiUAi 


12=1 1 

El 

2 




j ^ ^ ^ ^01 ■ 


(429) 
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The tunneling matrix elements Ja already discussed in section |2.4| as well as the atomic interaction param¬ 
eters 


Ua = 

J igs+9a){\W+i\^ + \W.i\^y -4ga\W+l\^\W.l\^ dx, 


Uo = 

Us= gs J \Woif dx , 

(430) 

K = 

{9s + 9a) J IWgil' (\W+lf + |W_,|") dx , 


C/a = 

^9aJ (Wo*,)"W+lW_ldx, 



and the relative shift of the mean energies of the eigenmodes 


6 = 


Ld 

kGlBZ 


4^)(k)-4°)(k) 


(431) 


can be simultaneously changed by varying the laser intensity and/or the angle 9, but the variations of 
and 5 are much faster. The parameter 17^ can be either positive or negative depending on the sign of the 
antisymmetric coupling ga- We will consider the case of repulsive interactions when [/g and C/a are positive 
and of about equal size but K can be larger or smaller than depending on the sign of the antisymmetric 
coupling ga- 


12.5.1. "Ferromagnetic" and "antiferromagnetic" superfluid states 


In the case when only A-mode is populated the Hamiltonian (429 \ becomes equivalent to that of spinless 
bosons. The only difference is that the tunneling matrix element J\ can take not only positive but also 
negative values. This leads to two different superfluid regimes which can be easily understood in the limit 
of the ideal gas. In the usual situation of positive Ja the eigenstates of non-interacting bosons are given 
by Eqs. (1491, (1501 and the ground state corresponds to k = 0, i.e., the phases of the wavefunction for 
any lattice site are the same. In the case of negative Ja, the sign in Eq. (1491 is reversed but Eq. (1501 
remains unchanged. Now the ground state corresponds to = tt which leads to the phase shift of the 
wavefunction for the neighboring lattice sites. In analogy to spin ordering in magnetic systems, one can call 
this “ferromagnetic" and “antiferromagnetic" phase ordering |489| . 

In the presence of nonvanishing interactions this qualitative picture is preserved and the two different 
superfluid phases are readily distinguishable experimentally via the spatial interference pattern generated by 
the coherent matter waves which one obtains in the time-of-flight images: The interference maxima obtained 
in the ferromagnetic case turn into minima in the antiferromagnetic case and vice versa. Deep in the MI 
phase, the phase coherence is lost and the sign of Ja does not play any role. 

The phase diagram of the system is determined by the ratio Ja/Ua. In the case of spinless bosons 
the ratio of the tunneling rate to the interaction parameter is a monotonic function of |Vo|- In the two- 
component case with the A-coupling we are dealing with, it has quite different properties. Its dependence on 
the control parameters Vg and 9 is not monotonic and there is a change of sign. Therefore, it is reasonable 
to draw g — Vq — 9 diagrams instead of /i — J diagrams. At the points (Vg, 9), where Ja/Ua vanishes, we 
have the MI phase for any values of g/U. These points define the boundary between the ferromagnetic and 
antiferromagnetic states. The mean-field phase diagram in the plane spanned by 9 and Vg for (/i/C/)c = ^2—1 
corresponding to the tip of the MI lobe for (hi) = 1 is shown in Fig. 61 The MI phase is strongly suppressed 


in the case Vg > 0 due to the dominant contribution of the dark component into the ground state. 


12.5.2. First- and second-order phase transitions 

In the present section, we are interested in the situations when both modes are occupied. Therefore we 
have to restrict ourselves to negative Vg and small values of 9. In this case, Ja as well as Jg are positive 
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Figure 61: Phase diagram in the {6, Vq) plane for {fi/U)c = \/2 — 1. Vq < 0: The boundaries between the SF and the MI 
({h\) = 1) phases are shown by solid lines. The dashed line corresponds to Ja = 0. It lies always in the MI phase, separating 
the regions of the ferromagnetic (SFq) and antiferromagnetic (SFtt) superfluid phases. Vq > 0: The line Ja = 0 as well as the 
two boundaries separating SFq and SF^^ superfiuid phases from the Mott phase are indistinguishable on the large-scale plot. 
The MI phase is located in the extremely narrow region between the SF phases. 


quantities. If ^ = 0, the A and 0 modes are degenerate and S vanishes. In addition, the components of the 
Wannier spinors satisfy the relation Hfyi = W-\ = VFoi/v^ and we have 

Jo = Ja = J,Uo = Ua = Us,K = U, + U^, (432) 


i.e., the Hamiltonian becomes symmetric with respect to the exchange of the indices 0 and A of the bosonic 
operators. Very useful representation of the Hamiltonian can be obtained in this case with the aid of the 
isospin operator Si with the components 


which has the property 


where 


<^11 = 2 (®A1®01 + ®01®Al) ’ 

S 21 = 2 (®A1®01 ~ ®01®Al) ’ 

*^31 = 2 (®01®01 ~ ®Al®Al) > 


S? = ^ ^ + 1 


fli — OoiOq, 


+ «Ah 


Al 


(433) 


(434) 

(435) 


is an operator of the total particle number on site 1. Note that the components of the operator Si generate 
the SU(2) algebra. In this notations, the Hamiltonian takes the form 


.Hbh = —J 


EEE( 

a iy—1 1 


at, a 


crl“CTl+e„ 


u. 




+ h.c. j + — ^ hi(ni — 1) — ^ hi + 2C7a ^ <5ii. 


(436) 


If the tunneling is negligible, the local eigenstates of the Hamiltonian (4361 coincide with the eigenstates 
|n/2, Al) of the isospin operator with the corresponding eigenenergies given by 


</2.A^ = - 1) + C/a (2M^ - ”) , (437) 

where n is the number of atoms on a lattice site and A4 = —nj2 ^... ,n/2 is the isospin projection on the 
direction 1 in the isospin space. Note that n/2 plays the role of the isospin quantum number. If n is even, 
the states with A4 = 0 are unique, while the others with Al yf 0 are doubly degenerate. If n is odd, the 
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Figure 62: (color online) Free energy of the Hamiltonian l |436^ for {Un/Us ~ —0.009 14581 '!. fi/Ua = 2.5. 2dJjUa = 

0.03 (a), 0.09 (b). 


states with = 0 do not exist and, therefore, all the eigenstates in the absence of tunneling are doubly 
degenerate. 

If C4 < 0, the ground state is described by = ±n/2. In the case of 14 > 0, one has to distinguish 
between odd and even n, where the ground state corresponds to = 0 and Ai = ±1/2, respectively. In 
the basis of the eigenstates of the isospin operator, nonvanishing matrix elements of the bosonic operators 
are given by 


- 1 .. 1 


n . A 

In - 1 ,, 1 


2 ’-^^2 

OA 

2 '-^) 

- 2 ’^^2 

oq 




■\/n ± 2 Ai 


■exp(T*0 


(438) 


Zero-temperature mean-field phase diagram of the symmetric spin-1/2 Bose-Hubbard model described 
by the Hamiltonian (4361 was studied in details in Refs. [1261149(114492] . The formalism is based on the 


decoupling approximation similar to Eq. (318) 


trli C7L2 




(439) 


where ^o-i = (oai) is the order parameter for Bose-Einstein condensation in the component cr = 0, A, which 
can be considered as a real and position-independent quantity. The free energy per lattice site A is in 
general independent of the sign of ipcr, ±'('0Aj V'o) = -±(|(/a|, 1401)> it is a symmetric function of 4o and 

4a: 7^(40, 4a) = A(4a,4o)- 

Typical dependences of 7 on the order parameters 4a and 4o in the case [/a < 0 are shown in Fig. 

As long as the ratio J/C4 is small, there is a single minimum at 4o = 4 a = 0 corresponding to the MI 
phase. For larger values of J/Ug the single minimum transforms into four equal minima located on the lines 
4a = ±4o- This corresponds to the SF phase with equal contributions of both components. As in the case 
of spinless bosons, the transition is continuous, i.e., second order, and the phase boundary is described by 


Eq. (329) with U = Us- |t7a 


The case of positive t/a is quite different as one can see in Figs. which shows typical dependences of 
the ground-state energy per lattice site on the order parameters. At small values of J/Ug, there is again only 
one minimum at 4a = 4o = 0 which corresponds to the MI phase. If the ratio J/Ug is increased, four equal 
minima appear on the lines 4a = 0 and 4o = 0 which means that the superfluid is polarized. In addition, 
the minimum at 4a = 4o = 0 does not always disappear immediately, but only if J/Ug is further increased. 
This implies that in such cases the phase transition is discontinuous, i.e., of the first order, and in a certain 
range of J/Ug the two phases coexist. 
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Figure 63; (color online) Free energy of the Hamiltonian ( |436| l for ^^Na {Ua/Us ~ 0.037 |458p . (a) fi/Us = 2.5, 2dJ/Us = 0.12. 
(b) V’A = 0, fi/Us = 1.5, 2dJ/Us = 0.125(1), 0.148(2), 0.157(^0.167(4). 
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Figure 64: In the shaded regions of this diagram .F(i/’) has two minima at certain values of J. In the remaining part .F('0) 
has only one minimum. (Reproduced from Ref. liSOl . ©2005 American Physical Society.) 


The values of and [ 4 , which allow the first-order transition are shown in Fig. ([64|. For n = 1, the 
transition is always second order. First-order transition is possible for n > 2 and U^/Us smaller than some 
critical value, which is about 0.188 for even n and grows from 0.012 (n = 3) to 0.015 (n — j. 00) for odd n. 
In the case of ^^Na shown in Fig. 65 an interesting regime is achieved, when the QPT for odd n is second 
order, but for even n it is first order. 


The phase diagram for ^^Na is presented in Fig. 65 It consists of a series of (internal) lobes corresponding 


to the stable Mott phase and external regions corresponding to the stable superfiuid phase. However, in the 
case of even n, the two regions are separated from one another by intermediate ones, where the stable and 
metastable superfiuid and Mott phases coexist. The boundary separating the region of the stable superfiuid 
phase from other ones can be calculated exactly in the mean-field theory using the second-order perturbation 
theory. In the case of odd n, the boundary is given by 


f = 4 


n — 1 


+ 


n + 3 


fi' — n + 1 + n + U^ — 


— -|- 2(?r -f 1) 


1 


n — U'^ — fi' fi' — n + 1J _ 


(440) 


where J' = 2dJ/Us, U'^ 
by the equation 


Ua/Us, and = ^/Us with Us{n — 1) < fi < UsTl — 14- For even n, it is described 


{n' -n+l + U'J {n - n') 

f^-n + l + U!, + (l + m^n/2 ' 


(441) 


where Us{n — 1) — U^ < H < C4n. 

The predictions of the mean-field theory have been tested by QMC calculations with one and two bosons 
per site in one and two dimensions |4911149,3] . It was verified that for ferromagnetic interactions (14 < 0) 
the phases are unpolarized and the transitions are continuous. In the case of antiferromagnetic interactions 
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Figure 65: (color online) Phase diagram for ^^Na {Ua/Us ~ 0.037 058]). The regions of metastable SF phase coexisting with 
the stable MI phase and that of metastable MI phase coexisting with the stable SF phase are denoted by MSF and MMI, 
respectively. (Adapted from Ref. liool . ©2005 American Physical Society.) 


(C/a > 0) it was confirmed that the SF phase is always dominated by one component. The first-order 
transitions were found in two dimensions with two bosons per lattice site and for t/a/C/g ^ 0.25 which is in 
good agreement with the results shown in Fig. However, the mean-field theory fails to describe all the 
features of the MI phase at finite tunneling and of the one-dimensional systems in general. QMC calculations 
demonstrated that all the transitions in one dimension are continuous. In addition, it was shown that the 
populations of the components in the case of one atoms per site become unbalanced already in the MI 
phase at finite ratios J/Us and this polarization persists through the MI-SF transition. It was also shown 
that thermal fluctuations immediately destroy the polarization of the MI phase, while the SF phase is less 
sensitive to that at least for small temperatures IMS]. 


13. Concluding remarks 

In this review, we presented studies of equilibrium properties of ultracold bosons with short-range in¬ 
teractions in optical lattices of the simplest hypercubic geometries. Although this kind of systems displays 
already quite rich physics, this is just a tiny part of activities in a huge area of research on cold atoms in 
optical lattices. Other topics include atoms with long-range dipole-dipole interactions, disordered systems, 
mixtures of different bosonic species, fermions as well as Bose-Fermi mixtures, lattices with more complicated 
geometric structures, e.g., hexagonal and triangular. We just mention intriguing experimental observations 
of Anderson localization and Bose-glass phase in incommensurate lattices |494[ 1495] . observation of effective 
multi-body interactions up to the six-body case in a three-dimensional lattice IMS], experimental realization 
of strong effective magnetic fields in a two-dimensional optical super lattice |497l 1498] , quantum simulations 
of a one-dimensional chain of interacting Ising spins in the presence of longitudinal and transverse fields IMS], 
in situ studies of photon-assisted tunneling in a Mott insulator by amplitude modulation of a tilted opti¬ 
cal lattice IMS], realization of a finite-momentum superfluid in the lowest P-band of a square lattice with 
two different depths of the potential wells arranged in a checkerboard pattern |501| . studies of quantum 
magnetism in high-spin systems |33| . observation of spin-exchange interactions with polar molecules |M^ . 

Another field of research which became very popular in the last years is the study of nonequilibrium 
phenomena. This is due to the fact that ultracold atoms are very well isolated from the environment which 
makes possible experimental investigations of relaxation and thermalization in closed quantum systems |503F 
IMS] as well as the dynamics of nonlocal quantum correlations |195| . This is a challenge for theorists because 
exact calculations of the dynamics of interacting quantum systems remain a difficult problem. In spite of 
a great progress achieved for one-dimensional systems with the aid of matrix-product states, the methods 
for higher-dimensional systems are not so well developed. Some progress in this direction is achieved by the 
dynamical mean-field theory |M1 150614508] and Monte Carlo methods |50911510] . An alternative approach 
was recently suggested in Refs. |401L 140411405L I5TT] which deals with the dynamical equations for the reduced 
density matrices of different number of lattice sites. However, the full theoretical description of sufficiently 
large systems in the whole parameter range and for arbitrarily long times is still not reached. 
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